module cam_history
  !-------------------------------------------------------------------------------------------
  !
  ! The cam_history module provides the user interface for CAM's history output capabilities.
  ! It maintains the lists of fields that are written to each history file, and the associated
  ! metadata for those fields such as descriptive names, physical units, time axis properties,
  ! etc.  It also contains the programmer interface which provides routines that are called from
  ! the physics and dynamics initialization routines to define the fields that are produced by
  ! the model and are available for output, and the routine that is called from the corresponding
  ! run method to add the field values into a history buffer so that they may be output to disk.
  !
  ! There are two special history files.  The initial file and the satellite track file.
  !
  ! Public functions/subroutines:
  !   addfld, add_default
  !   intht
  !   history_initialized
  !   write_restart_history
  !   read_restart_history
  !   outfld
  !   wshist
  !-----------------------------------------------------------------------

   use shr_kind_mod,    only: r8 => shr_kind_r8, r4 => shr_kind_r4
   use shr_kind_mod,    only: cl=>SHR_KIND_CL
   use shr_sys_mod,     only: shr_sys_flush
   use spmd_utils,      only: masterproc
   use ppgrid,          only: pcols, psubcols
   use cam_instance,    only: inst_suffix
   use cam_control_mod, only: caseid, ctitle
   use filenames,       only: interpret_filename_spec
   use cam_initfiles,   only: ncdata, bnd_topo
   use cam_abortutils,  only: endrun

   use pio,          only: file_desc_t, var_desc_t, pio_setframe, pio_write,  &
                           pio_noerr, pio_bcast_error, pio_internal_error,    &
                           pio_seterrorhandling, pio_get_var, pio_clobber,    &
                           pio_int, pio_real, pio_double, pio_char,           &
                           pio_offset_kind, pio_unlimited, pio_global,        &
                           pio_inq_dimlen, pio_def_var, pio_enddef,           &
                           pio_put_att, pio_put_var, pio_get_att


   use perf_mod,            only: t_startf, t_stopf
   use cam_logfile,         only: iulog
   use cam_history_support, only: max_fieldname_len, fieldname_suffix_len,    &
                                  max_chars, ptapes, fieldname_len,           &
                                  max_string_len, date2yyyymmdd, pflds,       &
                                  fieldname_lenp2, sec2hms,                   &
                                  field_info, active_entry, hentry,           &
                                  horiz_only, write_hist_coord_attrs,         &
                                  write_hist_coord_vars, interp_info_t,       &
                                  lookup_hist_coord_indices, get_hist_coord_index
   use sat_hist,            only: is_satfile
   use solar_parms_data,    only: solar_parms_on, kp=>solar_parms_kp, ap=>solar_parms_ap
   use solar_parms_data,    only: f107=>solar_parms_f107, f107a=>solar_parms_f107a, f107p=>solar_parms_f107p
   use solar_wind_data,     only: solar_wind_on, byimf=>solar_wind_byimf, bzimf=>solar_wind_bzimf
   use solar_wind_data,     only: swvel=>solar_wind_swvel, swden=>solar_wind_swden

  implicit none
  private
  save

  ! Forward common parameters to present unified interface to cam_history
  public :: fieldname_len, horiz_only
  public :: get_field_properties
  public :: cam_history_snapshot_deactivate
  public :: cam_history_snapshot_activate

  !
  ! master_entry: elements of an entry in the master field list
  !
  type master_entry
    type (field_info)        :: field            ! field information
    character(len=max_fieldname_len) :: meridional_field = '' ! for vector fields
    character(len=max_fieldname_len) :: zonal_field = '' ! for vector fields
    character(len=1)         :: avgflag(ptapes)  ! averaging flag
    character(len=max_chars) :: time_op(ptapes)  ! time operator (e.g. max, min, avg)
    logical               :: act_sometape     ! Field is active on some tape
    logical               :: actflag(ptapes)  ! Per tape active/inactive flag
    integer               :: htapeindx(ptapes)! This field's index on particular history tape
    type(master_entry), pointer :: next_entry => null() ! The next master entry
  end type master_entry

  type (master_entry), pointer :: masterlinkedlist => null()   ! master field linkedlist top

  type master_list
    type(master_entry), pointer :: thisentry => null()
  end type master_list

  type (master_list), pointer :: masterlist(:) => null() ! master field array for hash lookup of field

  ! history tape info
  type (active_entry), pointer :: tape(:) => null()          ! history tapes
  type (active_entry), target,allocatable :: history_tape(:) ! history tapes
  type (active_entry), target, allocatable :: restarthistory_tape(:) ! restart history tapes

  type rvar_id
    type(var_desc_t), pointer :: vdesc => null()
    integer :: type
    integer :: ndims
    integer :: dims(4)
    character(len=fieldname_lenp2) :: name
  end type rvar_id
  type rdim_id
    integer :: len
    integer :: dimid
    character(len=fieldname_lenp2) :: name
  end type rdim_id
  !
  !   The size of these parameters should match the assignments in restart_vars_setnames and restart_dims_setnames below
  !
  integer, parameter :: restartvarcnt              = 38
  integer, parameter :: restartdimcnt              = 10
  type(rvar_id)      :: restartvars(restartvarcnt)
  type(rdim_id)      :: restartdims(restartdimcnt)
  integer, parameter :: ptapes_dim_ind             =  1
  integer, parameter :: max_string_len_dim_ind     =  2
  integer, parameter :: fieldname_lenp2_dim_ind    =  3
  integer, parameter :: pflds_dim_ind              =  4
  integer, parameter :: max_chars_dim_ind          =  5
  integer, parameter :: max_fieldname_len_dim_ind  =  6
  integer, parameter :: maxnflds_dim_ind           =  7
  integer, parameter :: maxvarmdims_dim_ind        =  8
  integer, parameter :: registeredmdims_dim_ind    =  9
  integer, parameter :: max_hcoordname_len_dim_ind = 10

  integer :: nfmaster = 0             ! number of fields in master field list
  integer :: nflds(ptapes)            ! number of fields per tape

  ! per tape sampling frequency (0=monthly avg)

  integer :: idx                      ! index for nhtfrq initialization
  integer :: nhtfrq(ptapes) = (/0, (-24, idx=2,ptapes)/)  ! history write frequency (0 = monthly)
  integer :: mfilt(ptapes) = 30       ! number of time samples per tape
  integer :: nfils(ptapes)            ! Array of no. of files on current h-file
  integer :: ndens(ptapes) = 2        ! packing density (double (1) or real (2))
  integer :: ncprec(ptapes) = -999    ! netcdf packing parameter based on ndens
  real(r8) :: beg_time(ptapes)        ! time at beginning of an averaging interval

  logical :: rgnht(ptapes) = .false.  ! flag array indicating regeneration volumes
  logical :: hstwr(ptapes) = .false.  ! Flag for history writes
  logical :: empty_htapes  = .false.  ! Namelist flag indicates no default history fields
  logical :: htapes_defined = .false. ! flag indicates history contents have been defined

  character(len=cl) :: model_doi_url = '' ! Model DOI
  ! NB: This name must match the group name in namelist_definition.xml
  character(len=*), parameter   :: history_namelist = 'cam_history_nl'
  character(len=max_string_len) :: hrestpath(ptapes) = (/(' ',idx=1,ptapes)/) ! Full history restart pathnames
  character(len=max_string_len) :: nfpath(ptapes) = (/(' ',idx=1,ptapes)/) ! Array of first pathnames, for header
  character(len=max_string_len) :: cpath(ptapes)                   ! Array of current pathnames
  character(len=max_string_len) :: nhfil(ptapes)                   ! Array of current file names
  character(len=1)  :: avgflag_pertape(ptapes) = (/(' ',idx=1,ptapes)/) ! per tape averaging flag
  character(len=16)  :: logname             ! user name
  character(len=16)  :: host                ! host name
  character(len=8)   :: inithist = 'YEARLY' ! If set to '6-HOURLY, 'DAILY', 'MONTHLY' or
  ! 'YEARLY' then write IC file
  logical            :: inithist_all = .false. ! Flag to indicate set of fields to be
                                          ! included on IC file
                                          !  .false.  include only required fields
                                          !  .true.   include required *and* optional fields
  character(len=fieldname_lenp2) :: fincl(pflds,ptapes) ! List of fields to add to primary h-file
  character(len=max_chars)       :: fincllonlat(pflds,ptapes) ! List of fields to add to primary h-file
  character(len=fieldname_lenp2) :: fexcl(pflds,ptapes) ! List of fields to rm from primary h-file
  character(len=fieldname_lenp2) :: fwrtpr(pflds,ptapes) ! List of fields to change default history output prec
  character(len=fieldname_suffix_len ) :: fieldname_suffix = '&IC' ! Suffix appended to field names for IC file

  ! Parameters for interpolated output tapes
  logical, public       :: interpolate_output(ptapes) = .false.
  ! The last two history files are not supported for interpolation
  type(interp_info_t)   :: interpolate_info(ptapes - 2)

  ! Allowed history averaging flags
  ! This should match namelist_definition.xml => avgflag_pertape (+ ' ')
  ! The presence of 'ABI' and 'XML' in this string is a coincidence
  character(len=7), parameter    :: HIST_AVG_FLAGS = ' ABIXML'
  character(len=22) ,parameter   :: LT_DESC = 'mean (over local time)' ! local time description
  logical :: collect_column_output(ptapes)

  integer :: maxvarmdims=1
  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  !  Hashing.
  !
  !  Accelerate outfld processing by using a hash function of the field name
  !  to index masterlist and determine whehter the particular field is to
  !  be written to any history tape.
  !
  !
  !  Note: the outfld hashing logic will fail if any of the following are true:
  !
  !         1) The lower bound on the dimension of 'masterlist' is less than 1.
  !
  !         2) 'outfld' is called with field names that are not defined on
  !            masterlist.  This applies to both initial/branch and restart
  !            runs.
  !
  !         3) An inconsistency between a field's tape active flag
  !            'masterlist(ff)%actflag(t)' and active fields read from
  !            restart files.
  !
  !         4) Invoking function 'gen_hash_key' before the primary and secondary
  !            hash tables have been created (routine bld_outfld_hash_tbls).
  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  !
  !  User definable constants for hash and overflow tables.
  !  Define size of primary hash table (specified as 2**size).
  !
  integer, parameter :: tbl_hash_pri_sz_lg2 = 16
  !
  !  Define size of overflow hash table % of primary hash table.
  !
  integer, parameter :: tbl_hash_oflow_percent = 20
  !
  !  Do *not* modify the parameters below.
  !
  integer, parameter :: tbl_hash_pri_sz = 2**tbl_hash_pri_sz_lg2
  integer, parameter :: tbl_hash_oflow_sz = tbl_hash_pri_sz * (tbl_hash_oflow_percent/100.0_r8)
  !
  !  The primary and overflow tables are organized to mimimize space (read:
  !  try to maximimze cache line usage).
  !
  !  gen_hash_key(fieldname) will return an index on the interval
  !  [0 ... tbl_hash_pri_sz-1].
  !
  !
  !  Primary:
  !  gen_hash_key(fieldname)-------+     +----------+
  !                                |     |   -ii    | 1 ------>tbl_hash_oflow(ii)
  !                                |     +----------+
  !                                +-->  |    ff    | 2 ------>masterlist(ff)
  !                                      +----------+
  !                                      |          | ...
  !                                      +----------+
  !                                      |          | tbl_hash_pri_sz
  !                                      +----------+
  !
  !  Overflow (if tbl_hash_pri() < 0):
  !  tbl_hash_pri(gen_hash_key(fieldname))
  !                         |
  !                         |            +----------+
  !                         |            |     1    | 1  (one entry on O.F. chain)
  !                         |            +----------+
  !                         |            |    ff_m  | 2
  !                         |            +----------+
  !                         +--------->  |     3    | 3  (three entries on chain)
  !                                      +----------+
  !                                      |    ff_x  | 4
  !                                      +----------+
  !                                      |    ff_y  | 5
  !                                      +----------+
  !                                      |    ff_z  | 6
  !                                      +----------+
  !                                      |          | ...
  !                                      +----------+
  !                                      |          | tbl_hash_oflow_sz
  !                                      +----------+
  !
  !
  integer, dimension(0:tbl_hash_pri_sz-1) :: tbl_hash_pri ! Primary hash table
  integer, dimension(tbl_hash_oflow_sz) :: tbl_hash_oflow ! Overflow hash table
  !
  !  Constants used in hashing function gen_hash_key.
  !  Note: if the constants in table 'tbl_gen_hash_key' below are modified,
  !        changes are required to routine 'gen_hash_key' because of specific
  !        logic in the routine that optimizes character strings of length 8.
  !

  integer, parameter :: gen_hash_key_offset = z'000053db'

  integer, parameter :: tbl_max_idx = 15  ! 2**N - 1
  integer, dimension(0:tbl_max_idx) :: tbl_gen_hash_key = &
       (/61,59,53,47,43,41,37,31,29,23,17,13,11,7,3,1/)

  !
  ! Filename specifiers for history, initial files and restart history files
  ! (%c = caseid, $y = year, $m = month, $d = day, $s = seconds in day, %t = tape number)
  !
  character(len=max_string_len) :: rhfilename_spec = '%c.cam.rh%t.%y-%m-%d-%s.nc' ! history restart
  character(len=max_string_len) :: hfilename_spec(ptapes) = (/ (' ', idx=1, ptapes) /) ! filename specifyer


  interface addfld
    module procedure addfld_1d
    module procedure addfld_nd
  end interface

  ! Needed by cam_diagnostics
  public :: inithist_all

  integer :: lcltod_start(ptapes) ! start time of day for local time averaging (sec)
  integer :: lcltod_stop(ptapes)  ! stop time of day for local time averaging, stop > start is wrap around (sec)

  ! Needed by stepon and cam_restart
  public :: hstwr
  public :: nfils, mfilt

  ! Functions
  public :: history_readnl            ! Namelist reader for CAM history
  public :: init_restart_history      ! Write restart history data
  public :: write_restart_history     ! Write restart history data
  public :: read_restart_history      ! Read restart history data
  public :: wshist                    ! Write files out
  public :: outfld                    ! Output a field
  public :: intht                     ! Initialization
  public :: history_initialized       ! .true. iff cam history initialized
  public :: wrapup                    ! process history files at end of run
  public :: write_inithist            ! logical flag to allow dump of IC history buffer to IC file
  public :: addfld                    ! Add a field to history file
  public :: add_default               ! Add the default fields
  public :: register_vector_field     ! Register vector field set for interpolated output
  public :: get_hfilepath             ! Return history filename
  public :: get_ptapes                ! Return the number of tapes being used
  public :: get_hist_restart_filepath ! Return the full filepath to the history restart file
  public :: hist_fld_active           ! Determine if a field is active on any history file
  public :: hist_fld_col_active       ! Determine if a field is active on any history file at
  ! each column in a chunk

CONTAINS

  subroutine intht (model_doi_url_in)
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Initialize history file handler for initial or continuation run.
    !          For example, on an initial run, this routine initializes "ptapes"
    !          history files.  On a restart or regeneration  run, this routine
    !          only initializes history files declared beyond what existed on the
    !          previous run.  Files which already existed on the previous run have
    !          already been initialized (i.e. named and opened) in routine RESTRT.
    !
    ! Method: Loop over tapes and fields per tape setting appropriate variables and
    !         calling appropriate routines
    !
    ! Author: CCM Core Group
    !
    !-----------------------------------------------------------------------
    use shr_sys_mod,      only: shr_sys_getenv
    use time_manager,     only: get_prev_time, get_curr_time
    use cam_control_mod,  only: restart_run, branch_run
    use sat_hist,         only: sat_hist_init
    use spmd_utils,       only: mpicom, masterprocid, mpi_character
    !
    !-----------------------------------------------------------------------
    !
    ! Dummy argument
    !
    character(len=cl), intent(in) :: model_doi_url_in
    !
    ! Local workspace
    !
    integer :: t, f              ! tape, field indices
    integer :: begdim1           ! on-node dim1 start index
    integer :: enddim1           ! on-node dim1 end index
    integer :: begdim2           ! on-node dim2 start index
    integer :: enddim2           ! on-node dim2 end index
    integer :: begdim3           ! on-node chunk or lat start index
    integer :: enddim3           ! on-node chunk or lat end index
    integer :: day, sec          ! day and seconds from base date
    integer :: rcode             ! shr_sys_getenv return code
    type(master_entry), pointer :: listentry
    character(len=32) :: fldname ! temp variable used to produce a left justified field name
    ! in the formatted logfile output

    !
    ! Save the DOI
    !
    model_doi_url = trim(model_doi_url_in)

    !
    ! Print master field list
    !

    if (masterproc) then
      write(iulog,*) ' '
      write(iulog,*)' ******* MASTER FIELD LIST *******'
    end if
    listentry=>masterlinkedlist
    f=0
    do while(associated(listentry))
      f=f+1
      if(masterproc) then
        fldname = listentry%field%name
        write(iulog,9000) f, fldname, listentry%field%units, listentry%field%numlev, &
             listentry%avgflag(1), trim(listentry%field%long_name)
9000    format(i5, 1x, a32, 1x, a16, 1x, i4, 1x, a1, 2x, a)
      end if
      listentry=>listentry%next_entry
    end do
    nfmaster = f
    if(masterproc) write(iulog,*)'intht:nfmaster=',nfmaster

    !
    !  Now that masterlinkedlist is defined and we are performing a restart run
    !  (after all addfld calls), construct primary and secondary hashing tables.
    !
    if (restart_run) then
       call print_active_fldlst()
       call bld_outfld_hash_tbls()
       call bld_htapefld_indices()
       return
    end if
    !
    ! Get users logname and machine hostname
    !
    if ( masterproc )then
      logname = ' '
      call shr_sys_getenv ('LOGNAME',logname,rcode)
      host = ' '
      call shr_sys_getenv ('HOST',host,rcode)
    end if
    ! PIO requires netcdf attributes have consistant values on all tasks
    call mpi_bcast(logname, len(logname), mpi_character, masterprocid, mpicom, rcode)
    call mpi_bcast(host,    len(host),    mpi_character, masterprocid, mpicom, rcode)
    !
    ! Override averaging flag for all fields on a particular tape if namelist input so specifies
    !
    do t=1,ptapes
      if (avgflag_pertape(t) /= ' ') then
        call h_override (t)
      end if
    end do
    !
    ! Define field list information for all history files.
    !
    call fldlst ()
    !
    ! Loop over max. no. of history files permitted
    !
    if (branch_run) then
      call get_prev_time(day, sec)  ! elapased time since reference date
    else
      call get_curr_time(day, sec)  ! elapased time since reference date
    end if
    do t=1,ptapes
      nfils(t) = 0            ! no. of time samples in hist. file no. t

      ! Time at beginning of current averaging interval.

      beg_time(t) = day + sec/86400._r8
    end do

    !
    ! Initialize history variables
    !
    do t=1,ptapes
      do f=1,nflds(t)
        begdim1  = tape(t)%hlist(f)%field%begdim1
        enddim1  = tape(t)%hlist(f)%field%enddim1
        begdim2  = tape(t)%hlist(f)%field%begdim2
        enddim2  = tape(t)%hlist(f)%field%enddim2
        begdim3  = tape(t)%hlist(f)%field%begdim3
        enddim3  = tape(t)%hlist(f)%field%enddim3
        allocate(tape(t)%hlist(f)%hbuf(begdim1:enddim1,begdim2:enddim2,begdim3:enddim3))
        tape(t)%hlist(f)%hbuf = 0._r8
        if (tape(t)%hlist(f)%avgflag .eq. 'S') then ! allocate the variance buffer for standard dev
           allocate(tape(t)%hlist(f)%sbuf(begdim1:enddim1,begdim2:enddim2,begdim3:enddim3))
           tape(t)%hlist(f)%sbuf = 0._r8
        endif
        if(tape(t)%hlist(f)%field%flag_xyfill .or. (avgflag_pertape(t) .eq. 'L')) then
          allocate (tape(t)%hlist(f)%nacs(begdim1:enddim1,begdim3:enddim3))
        else
          allocate (tape(t)%hlist(f)%nacs(1,begdim3:enddim3))
        end if
        tape(t)%hlist(f)%nacs(:,:) = 0
        tape(t)%hlist(f)%field%meridional_complement = -1
        tape(t)%hlist(f)%field%zonal_complement = -1
      end do
    end do
    ! Setup vector pairs for unstructured grid interpolation
    call setup_interpolation_and_define_vector_complements()
    !  Initialize the sat following history subsystem
    call sat_hist_init()

    return
  end subroutine intht

  logical function history_initialized()
    history_initialized = associated(masterlist)
  end function history_initialized

  subroutine history_readnl(nlfile)

    use namelist_utils, only: find_group_name
    use units,          only: getunit, freeunit
    use spmd_utils,     only: masterproc, masterprocid, mpicom
    use spmd_utils,     only: mpi_integer, mpi_logical, mpi_character
    use shr_string_mod, only: shr_string_toUpper
    use time_manager,   only: get_step_size
    use sat_hist,       only: sat_hist_readnl

    ! Dummy argument
    character(len=*), intent(in)   :: nlfile  ! filepath of namelist input file

    !
    ! Local variables
    integer                        :: dtime   ! Step time in seconds
    integer                        :: unitn, ierr, f, t
    character(len=8)               :: ctemp      ! Temporary character string

    character(len=fieldname_lenp2) :: fincl1(pflds)
    character(len=fieldname_lenp2) :: fincl2(pflds)
    character(len=fieldname_lenp2) :: fincl3(pflds)
    character(len=fieldname_lenp2) :: fincl4(pflds)
    character(len=fieldname_lenp2) :: fincl5(pflds)
    character(len=fieldname_lenp2) :: fincl6(pflds)
    character(len=fieldname_lenp2) :: fincl7(pflds)
    character(len=fieldname_lenp2) :: fincl8(pflds)
    character(len=fieldname_lenp2) :: fincl9(pflds)
    character(len=fieldname_lenp2) :: fincl10(pflds)

    character(len=max_chars)       :: fincl1lonlat(pflds)
    character(len=max_chars)       :: fincl2lonlat(pflds)
    character(len=max_chars)       :: fincl3lonlat(pflds)
    character(len=max_chars)       :: fincl4lonlat(pflds)
    character(len=max_chars)       :: fincl5lonlat(pflds)
    character(len=max_chars)       :: fincl6lonlat(pflds)
    character(len=max_chars)       :: fincl7lonlat(pflds)
    character(len=max_chars)       :: fincl8lonlat(pflds)
    character(len=max_chars)       :: fincl9lonlat(pflds)
    character(len=max_chars)       :: fincl10lonlat(pflds)

    character(len=fieldname_len)   :: fexcl1(pflds)
    character(len=fieldname_len)   :: fexcl2(pflds)
    character(len=fieldname_len)   :: fexcl3(pflds)
    character(len=fieldname_len)   :: fexcl4(pflds)
    character(len=fieldname_len)   :: fexcl5(pflds)
    character(len=fieldname_len)   :: fexcl6(pflds)
    character(len=fieldname_len)   :: fexcl7(pflds)
    character(len=fieldname_len)   :: fexcl8(pflds)
    character(len=fieldname_len)   :: fexcl9(pflds)
    character(len=fieldname_len)   :: fexcl10(pflds)

    character(len=fieldname_lenp2) :: fwrtpr1(pflds)
    character(len=fieldname_lenp2) :: fwrtpr2(pflds)
    character(len=fieldname_lenp2) :: fwrtpr3(pflds)
    character(len=fieldname_lenp2) :: fwrtpr4(pflds)
    character(len=fieldname_lenp2) :: fwrtpr5(pflds)
    character(len=fieldname_lenp2) :: fwrtpr6(pflds)
    character(len=fieldname_lenp2) :: fwrtpr7(pflds)
    character(len=fieldname_lenp2) :: fwrtpr8(pflds)
    character(len=fieldname_lenp2) :: fwrtpr9(pflds)
    character(len=fieldname_lenp2) :: fwrtpr10(pflds)

    integer                        :: interpolate_nlat(size(interpolate_info))
    integer                        :: interpolate_nlon(size(interpolate_info))
    integer                        :: interpolate_gridtype(size(interpolate_info))
    integer                        :: interpolate_type(size(interpolate_info))

    ! History namelist items
    namelist /cam_history_nl/ ndens, nhtfrq, mfilt, inithist, inithist_all,    &
         avgflag_pertape, empty_htapes, lcltod_start, lcltod_stop,             &
         fincl1lonlat, fincl2lonlat, fincl3lonlat, fincl4lonlat, fincl5lonlat, &
         fincl6lonlat, fincl7lonlat, fincl8lonlat, fincl9lonlat,               &
         fincl10lonlat, collect_column_output, hfilename_spec,                 &
         fincl1,  fincl2,  fincl3,  fincl4,  fincl5,                           &
         fincl6,  fincl7,  fincl8,  fincl9,  fincl10,                          &
         fexcl1,  fexcl2,  fexcl3,  fexcl4,  fexcl5,                           &
         fexcl6,  fexcl7,  fexcl8,  fexcl9,  fexcl10,                          &
         fwrtpr1, fwrtpr2, fwrtpr3, fwrtpr4, fwrtpr5,                          &
         fwrtpr6, fwrtpr7, fwrtpr8, fwrtpr9, fwrtpr10,                         &
         interpolate_nlat, interpolate_nlon,                                   &
         interpolate_gridtype, interpolate_type, interpolate_output

    ! Set namelist defaults (these should match initial values if given)
    fincl(:,:)               = ' '
    fincllonlat(:,:)         = ' '
    fexcl(:,:)               = ' '
    fwrtpr(:,:)              = ' '
    collect_column_output(:) = .false.
    avgflag_pertape(:)       = ' '
    ndens                    = 2
    nhtfrq(1)                = 0
    nhtfrq(2:)               = -24
    mfilt                    = 30
    inithist                 = 'YEARLY'
    inithist_all             = .false.
    empty_htapes             = .false.
    lcltod_start(:)          = 0
    lcltod_stop(:)           = 0
    hfilename_spec(:)        = ' '
    interpolate_nlat(:)      = 0
    interpolate_nlon(:)      = 0
    interpolate_gridtype(:)  = 1
    interpolate_type(:)      = 1
    interpolate_output(:)    = .false.

    ! Initialize namelist 'temporary variables'
    do f = 1, pflds
      fincl1(f)        = ' '
      fincl2(f)        = ' '
      fincl3(f)        = ' '
      fincl4(f)        = ' '
      fincl5(f)        = ' '
      fincl6(f)        = ' '
      fincl7(f)        = ' '
      fincl8(f)        = ' '
      fincl9(f)        = ' '
      fincl10(f)       = ' '
      fincl1lonlat(f)  = ' '
      fincl2lonlat(f)  = ' '
      fincl3lonlat(f)  = ' '
      fincl4lonlat(f)  = ' '
      fincl5lonlat(f)  = ' '
      fincl6lonlat(f)  = ' '
      fincl7lonlat(f)  = ' '
      fincl8lonlat(f)  = ' '
      fincl9lonlat(f)  = ' '
      fincl10lonlat(f) = ' '
      fexcl1(f)        = ' '
      fexcl2(f)        = ' '
      fexcl3(f)        = ' '
      fexcl4(f)        = ' '
      fexcl5(f)        = ' '
      fexcl6(f)        = ' '
      fexcl7(f)        = ' '
      fexcl8(f)        = ' '
      fexcl9(f)        = ' '
      fexcl10(f)       = ' '
      fwrtpr1(f)       = ' '
      fwrtpr2(f)       = ' '
      fwrtpr3(f)       = ' '
      fwrtpr4(f)       = ' '
      fwrtpr5(f)       = ' '
      fwrtpr6(f)       = ' '
      fwrtpr7(f)       = ' '
      fwrtpr8(f)       = ' '
      fwrtpr9(f)       = ' '
      fwrtpr10(f)      = ' '
    end do

    if (trim(history_namelist) /= 'cam_history_nl') then
      call endrun('HISTORY_READNL: CAM history namelist mismatch')
    end if
    if (masterproc) then
      write(iulog, *) 'Read in ',history_namelist,' namelist from: ',trim(nlfile)
      unitn = getunit()
      open(unitn, file=trim(nlfile), status='old')
      call find_group_name(unitn, history_namelist, status=ierr)
      if (ierr == 0) then
        read(unitn, cam_history_nl, iostat=ierr)
        if (ierr /= 0) then
          call endrun('history_readnl: ERROR reading namelist, '//trim(history_namelist))
        end if
      end if
      close(unitn)
      call freeunit(unitn)

      do f = 1, pflds
        fincl(f, 1) = fincl1(f)
        fincl(f, 2) = fincl2(f)
        fincl(f, 3) = fincl3(f)
        fincl(f, 4) = fincl4(f)
        fincl(f, 5) = fincl5(f)
        fincl(f, 6) = fincl6(f)
        fincl(f, 7) = fincl7(f)
        fincl(f, 8) = fincl8(f)
        fincl(f, 9) = fincl9(f)
        fincl(f,10) = fincl10(f)

        fincllonlat(f, 1) = fincl1lonlat(f)
        fincllonlat(f, 2) = fincl2lonlat(f)
        fincllonlat(f, 3) = fincl3lonlat(f)
        fincllonlat(f, 4) = fincl4lonlat(f)
        fincllonlat(f, 5) = fincl5lonlat(f)
        fincllonlat(f, 6) = fincl6lonlat(f)
        fincllonlat(f, 7) = fincl7lonlat(f)
        fincllonlat(f, 8) = fincl8lonlat(f)
        fincllonlat(f, 9) = fincl9lonlat(f)
        fincllonlat(f,10) = fincl10lonlat(f)

        fexcl(f, 1) = fexcl1(f)
        fexcl(f, 2) = fexcl2(f)
        fexcl(f, 3) = fexcl3(f)
        fexcl(f, 4) = fexcl4(f)
        fexcl(f, 5) = fexcl5(f)
        fexcl(f, 6) = fexcl6(f)
        fexcl(f, 7) = fexcl7(f)
        fexcl(f, 8) = fexcl8(f)
        fexcl(f, 9) = fexcl9(f)
        fexcl(f,10) = fexcl10(f)

        fwrtpr(f, 1) = fwrtpr1(f)
        fwrtpr(f, 2) = fwrtpr2(f)
        fwrtpr(f, 3) = fwrtpr3(f)
        fwrtpr(f, 4) = fwrtpr4(f)
        fwrtpr(f, 5) = fwrtpr5(f)
        fwrtpr(f, 6) = fwrtpr6(f)
        fwrtpr(f, 7) = fwrtpr7(f)
        fwrtpr(f, 8) = fwrtpr8(f)
        fwrtpr(f, 9) = fwrtpr9(f)
        fwrtpr(f,10) = fwrtpr10(f)
      end do

      !
      ! If generate an initial conditions history file as an auxillary tape:
      !
      ctemp = shr_string_toUpper(inithist)
      inithist = trim(ctemp)
      if ( (inithist /= '6-HOURLY') .and. (inithist /= 'DAILY')  .and.        &
           (inithist /= 'MONTHLY')  .and. (inithist /= 'YEARLY') .and.        &
           (inithist /= 'CAMIOP')   .and. (inithist /= 'ENDOFRUN')) then
        inithist = 'NONE'
      end if
      !
      ! History file write times
      ! Convert write freq. of hist files from hours to timesteps if necessary.
      !
      dtime = get_step_size()
      do t = 1, ptapes
        if (nhtfrq(t) < 0) then
          nhtfrq(t) = nint((-nhtfrq(t) * 3600._r8) / dtime)
        end if
      end do
      !
      ! Initialize the filename specifier if not already set
      ! This is the format for the history filenames:
      ! %c= caseid, %t=tape no., %y=year, %m=month, %d=day, %s=second, %%=%
      ! See the filenames module for more information
      !
      do t = 1, ptapes
        if ( len_trim(hfilename_spec(t)) == 0 )then
          if ( nhtfrq(t) == 0 )then
            ! Monthly files
            hfilename_spec(t) = '%c.cam' // trim(inst_suffix) // '.h%t.%y-%m.nc'
          else
            hfilename_spec(t) = '%c.cam' // trim(inst_suffix) // '.h%t.%y-%m-%d-%s.nc'
          end if
        end if
        !
        ! Only one time sample allowed per monthly average file
        !
        if (nhtfrq(t) == 0) then
          mfilt(t) = 1
        end if
      end do
    end if ! masterproc

    ! Print per-tape averaging flags
    if (masterproc) then
      do t = 1, ptapes
        if (avgflag_pertape(t) /= ' ') then
          write(iulog,*)'Unless overridden by namelist input on a per-field basis (FINCL),'
          write(iulog,*)'All fields on history file ',t,' will have averaging flag ',avgflag_pertape(t)
        end if
        ! Enforce no interpolation for satellite files
        if (is_satfile(t) .and. interpolate_output(t)) then
          write(iulog, *) 'WARNING: Interpolated output not supported for a satellite history file, ignored'
          interpolate_output(t) = .false.
        end if
        ! Enforce no interpolation for IC files
        if (is_initfile(t) .and. interpolate_output(t)) then
          write(iulog, *) 'WARNING: Interpolated output not supported for an initial data (IC) history file, ignored'
          interpolate_output(t) = .false.
        end if
      end do
    end if

    ! Write out inithist info
    if (masterproc) then
      if (inithist == '6-HOURLY' ) then
        write(iulog,*)'Initial conditions history files will be written 6-hourly.'
      else if (inithist == 'DAILY' ) then
        write(iulog,*)'Initial conditions history files will be written daily.'
      else if (inithist == 'MONTHLY' ) then
        write(iulog,*)'Initial conditions history files will be written monthly.'
      else if (inithist == 'YEARLY' ) then
        write(iulog,*)'Initial conditions history files will be written yearly.'
      else if (inithist == 'CAMIOP' ) then
        write(iulog,*)'Initial conditions history files will be written for IOP.'
      else if (inithist == 'ENDOFRUN' ) then
        write(iulog,*)'Initial conditions history files will be written at end of run.'
      else
        write(iulog,*)'Initial conditions history files will not be created'
      end if
    end if

    ! Print out column-output information
    do t = 1, size(fincllonlat, 2)
      if (ANY(len_trim(fincllonlat(:,t)) > 0)) then
        if (collect_column_output(t)) then
          write(iulog, '(a,i2,a)') 'History file, ',t,', has patch output, columns will be collected into ncol dimension'
        else
          write(iulog, '(a,i2,a)') 'History file, ',t,', has patch output, patches will be written to individual variables'
        end if
      end if
    end do

    ! Broadcast namelist variables
    call mpi_bcast(ndens, ptapes, mpi_integer, masterprocid, mpicom, ierr)
    call mpi_bcast(nhtfrq, ptapes, mpi_integer, masterprocid, mpicom, ierr)
    call mpi_bcast(mfilt, ptapes, mpi_integer, masterprocid, mpicom, ierr)
    call mpi_bcast(inithist,len(inithist), mpi_character, masterprocid, mpicom, ierr)
    call mpi_bcast(inithist_all,1, mpi_logical, masterprocid, mpicom, ierr)
    call mpi_bcast(lcltod_start, ptapes, mpi_integer, masterprocid, mpicom, ierr)
    call mpi_bcast(lcltod_stop,  ptapes, mpi_integer, masterprocid, mpicom, ierr)
    call mpi_bcast(collect_column_output, ptapes, mpi_logical, masterprocid, mpicom, ierr)
    call mpi_bcast(empty_htapes,1, mpi_logical, masterprocid, mpicom, ierr)
    call mpi_bcast(avgflag_pertape, ptapes, mpi_character, masterprocid, mpicom, ierr)
    call mpi_bcast(hfilename_spec, len(hfilename_spec(1))*ptapes, mpi_character, masterprocid, mpicom, ierr)
    call mpi_bcast(fincl, len(fincl (1,1))*pflds*ptapes, mpi_character, masterprocid, mpicom, ierr)
    call mpi_bcast(fexcl, len(fexcl (1,1))*pflds*ptapes, mpi_character, masterprocid, mpicom, ierr)

    call mpi_bcast(fincllonlat, len(fincllonlat (1,1))*pflds*ptapes, mpi_character, masterprocid, mpicom, ierr)

    call mpi_bcast(fwrtpr, len(fwrtpr(1,1))*pflds*ptapes, mpi_character, masterprocid, mpicom, ierr)
    t = size(interpolate_nlat, 1)
    call mpi_bcast(interpolate_nlat, t, mpi_integer, masterprocid, mpicom, ierr)
    call mpi_bcast(interpolate_nlon, t, mpi_integer, masterprocid, mpicom, ierr)
    call mpi_bcast(interpolate_gridtype, t, mpi_integer, masterprocid, mpicom, ierr)
    call mpi_bcast(interpolate_type, t, mpi_integer, masterprocid, mpicom, ierr)
    call mpi_bcast(interpolate_output, ptapes, mpi_logical, masterprocid, mpicom, ierr)

    ! Setup the interpolate_info structures
    do t = 1, size(interpolate_info)
      interpolate_info(t)%interp_type = interpolate_type(t)
      interpolate_info(t)%interp_gridtype = interpolate_gridtype(t)
      interpolate_info(t)%interp_nlat = interpolate_nlat(t)
      interpolate_info(t)%interp_nlon = interpolate_nlon(t)
    end do

    ! separate namelist reader for the satellite history file
    call sat_hist_readnl(nlfile, hfilename_spec, mfilt, fincl, nhtfrq, avgflag_pertape)

  end subroutine history_readnl

!==================================================================================================

  subroutine set_field_dimensions(field)
    use cam_history_support, only: hist_coord_size
    use cam_grid_support,    only: cam_grid_get_array_bounds, cam_grid_is_block_indexed
    ! Dummy arguments
    type(field_info), intent(inout) :: field

    ! Local variables
    integer                         :: i
    integer                         :: msize
    integer                         :: dimbounds(2,2)

    call cam_grid_get_array_bounds(field%decomp_type, dimbounds)
    field%begdim1  = dimbounds(1,1)
    field%enddim1  = dimbounds(1,2)
    field%begdim2  = 1
    if (associated(field%mdims)) then
      if (size(field%mdims) > 0) then
        field%enddim2  = 1
        do i = 1, size(field%mdims)
          msize = hist_coord_size(field%mdims(i))
          if (msize <= 0) then
            call endrun('set_field_dimensions: mdim size must be > 0')
          end if
          field%enddim2 = field%enddim2 * msize
        end do
      else
        if (field%numlev < 1) then
          if (masterproc) then
            write(iulog, *) 'SET_FIELD_DIMENSIONS WARNING: illegal numlev for ', trim(field%name)
          end if
          field%numlev = 1
        end if
        field%enddim2 = field%numlev
      end if
    else
      if (field%numlev < 1) then
        if (masterproc) then
          write(iulog, *) 'SET_FIELD_DIMENSIONS WARNING: illegal numlev for ', trim(field%name)
        end if
        field%numlev = 1
      end if
      field%enddim2 = field%numlev
    end if
    field%begdim3  = dimbounds(2,1)
    field%enddim3  = dimbounds(2,2)
    field%colperchunk = cam_grid_is_block_indexed(field%decomp_type)

  end subroutine set_field_dimensions

  subroutine setup_interpolation_and_define_vector_complements()
    use cam_interp_mod, only: setup_history_interpolation

    ! Local variables
    integer :: hf, f, ff
    logical :: interp_ok
    character(len=max_fieldname_len) :: mname
    character(len=max_fieldname_len) :: zname
    character(len=*), parameter      :: subname='setup_interpolation_and_define_vector_complements'

    ! Do not interpolate IC history and sat hist files
    if (any(interpolate_output)) then
      call setup_history_interpolation(interp_ok, ptapes-2,                   &
           interpolate_output, interpolate_info)
      do hf = 1, ptapes - 2
        if((.not. is_satfile(hf)) .and. (.not. is_initfile(hf))) then
          do f = 1, nflds(hf)
            if (field_part_of_vector(trim(tape(hf)%hlist(f)%field%name),      &
                 mname, zname)) then
              if (len_trim(mname) > 0) then
                ! This field is a zonal part of a set, find the meridional partner
                do ff = 1, nflds(hf)
                  if (trim(mname) == trim(tape(hf)%hlist(ff)%field%name)) then
                    tape(hf)%hlist(f)%field%meridional_complement = ff
                    tape(hf)%hlist(ff)%field%zonal_complement     = f
                    exit
                  end if
                  if (ff == nflds(hf)) then
                    call endrun(trim(subname)//': No meridional match for '//trim(tape(hf)%hlist(f)%field%name))
                  end if
                end do
              else if (len_trim(zname) > 0) then
                ! This field is a meridional part of a set, find the zonal partner
                do ff = 1, nflds(hf)
                  if (trim(zname) == trim(tape(hf)%hlist(ff)%field%name)) then
                    tape(hf)%hlist(f)%field%zonal_complement       = ff
                    tape(hf)%hlist(ff)%field%meridional_complement = f
                    exit
                  end if
                  if (ff == nflds(hf)) then
                    call endrun(trim(subname)//': No zonal match for '//trim(tape(hf)%hlist(f)%field%name))
                  end if
                end do
              else
                call endrun(trim(subname)//': INTERNAL ERROR, bad vector field')
              end if
            end if
          end do
        end if
      end do
    end if
  end subroutine setup_interpolation_and_define_vector_complements

  subroutine restart_vars_setnames()

    ! Local variable
    integer :: rvindex

    rvindex = 1
    restartvars(rvindex)%name = 'rgnht'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 1
    restartvars(rvindex)%dims(1) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'nhtfrq'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 1
    restartvars(rvindex)%dims(1) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'nflds'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 1
    restartvars(rvindex)%dims(1) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'nfils'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 1
    restartvars(rvindex)%dims(1) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'mfilt'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 1
    restartvars(rvindex)%dims(1) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'nfpath'
    restartvars(rvindex)%type = pio_char
    restartvars(rvindex)%ndims = 2
    restartvars(rvindex)%dims(1) = max_string_len_dim_ind
    restartvars(rvindex)%dims(2) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'cpath'
    restartvars(rvindex)%type = pio_char
    restartvars(rvindex)%ndims = 2
    restartvars(rvindex)%dims(1) = max_string_len_dim_ind
    restartvars(rvindex)%dims(2) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'nhfil'
    restartvars(rvindex)%type = pio_char
    restartvars(rvindex)%ndims = 2
    restartvars(rvindex)%dims(1) = max_string_len_dim_ind
    restartvars(rvindex)%dims(2) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'ndens'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 1
    restartvars(rvindex)%dims(1) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'fincllonlat'
    restartvars(rvindex)%type = pio_char
    restartvars(rvindex)%ndims = 3
    restartvars(rvindex)%dims(1) = max_chars_dim_ind
    restartvars(rvindex)%dims(2) = pflds_dim_ind
    restartvars(rvindex)%dims(3) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'ncprec'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 1
    restartvars(rvindex)%dims(1) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'beg_time'
    restartvars(rvindex)%type = pio_double
    restartvars(rvindex)%ndims = 1
    restartvars(rvindex)%dims(1) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'fincl'
    restartvars(rvindex)%type = pio_char
    restartvars(rvindex)%ndims = 3
    restartvars(rvindex)%dims(1) = fieldname_lenp2_dim_ind
    restartvars(rvindex)%dims(2) = pflds_dim_ind
    restartvars(rvindex)%dims(3) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'fexcl'
    restartvars(rvindex)%type = pio_char
    restartvars(rvindex)%ndims = 3
    restartvars(rvindex)%dims(1) = fieldname_lenp2_dim_ind
    restartvars(rvindex)%dims(2) = pflds_dim_ind
    restartvars(rvindex)%dims(3) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'field_name'
    restartvars(rvindex)%type = pio_char
    restartvars(rvindex)%ndims = 3
    restartvars(rvindex)%dims(1) = max_fieldname_len_dim_ind
    restartvars(rvindex)%dims(2) = maxnflds_dim_ind
    restartvars(rvindex)%dims(3) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'decomp_type'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 2
    restartvars(rvindex)%dims(1) = maxnflds_dim_ind
    restartvars(rvindex)%dims(2) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'numlev'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 2
    restartvars(rvindex)%dims(1) = maxnflds_dim_ind
    restartvars(rvindex)%dims(2) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'hrestpath'
    restartvars(rvindex)%type = pio_char
    restartvars(rvindex)%ndims = 2
    restartvars(rvindex)%dims(1) = max_string_len_dim_ind
    restartvars(rvindex)%dims(2) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'hwrt_prec'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 2
    restartvars(rvindex)%dims(1) = maxnflds_dim_ind
    restartvars(rvindex)%dims(2) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'avgflag'
    restartvars(rvindex)%type = pio_char
    restartvars(rvindex)%ndims = 2
    restartvars(rvindex)%dims(1) = maxnflds_dim_ind
    restartvars(rvindex)%dims(2) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'sampling_seq'
    restartvars(rvindex)%type = pio_char
    restartvars(rvindex)%ndims = 3
    restartvars(rvindex)%dims(1) = max_chars_dim_ind
    restartvars(rvindex)%dims(2) = maxnflds_dim_ind
    restartvars(rvindex)%dims(3) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'cell_methods'
    restartvars(rvindex)%type = pio_char
    restartvars(rvindex)%ndims = 3
    restartvars(rvindex)%dims(1) = max_chars_dim_ind
    restartvars(rvindex)%dims(2) = maxnflds_dim_ind
    restartvars(rvindex)%dims(3) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'long_name'
    restartvars(rvindex)%type = pio_char
    restartvars(rvindex)%ndims = 3
    restartvars(rvindex)%dims(1) = max_chars_dim_ind
    restartvars(rvindex)%dims(2) = maxnflds_dim_ind
    restartvars(rvindex)%dims(3) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'units'
    restartvars(rvindex)%type = pio_char
    restartvars(rvindex)%ndims = 3
    restartvars(rvindex)%dims(1) = max_chars_dim_ind
    restartvars(rvindex)%dims(2) = maxnflds_dim_ind
    restartvars(rvindex)%dims(3) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'xyfill'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 2
    restartvars(rvindex)%dims(1) = maxnflds_dim_ind
    restartvars(rvindex)%dims(2) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'lcltod_start'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 1
    restartvars(rvindex)%dims(1) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'lcltod_stop'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 1
    restartvars(rvindex)%dims(1) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'fillvalue'
    restartvars(rvindex)%type = pio_double
    restartvars(rvindex)%ndims = 2
    restartvars(rvindex)%dims(1) = maxnflds_dim_ind
    restartvars(rvindex)%dims(2) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'mdims'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 3
    restartvars(rvindex)%dims(1) = maxvarmdims_dim_ind
    restartvars(rvindex)%dims(2) = maxnflds_dim_ind
    restartvars(rvindex)%dims(3) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'mdimnames'
    restartvars(rvindex)%type = pio_char
    restartvars(rvindex)%ndims = 2
    restartvars(rvindex)%dims(1) = max_hcoordname_len_dim_ind
    restartvars(rvindex)%dims(2) = registeredmdims_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'is_subcol'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 2
    restartvars(rvindex)%dims(1) = maxnflds_dim_ind
    restartvars(rvindex)%dims(2) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'interpolate_output'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 1
    restartvars(rvindex)%dims(1) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'interpolate_type'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 1
    restartvars(rvindex)%dims(1) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'interpolate_gridtype'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 1
    restartvars(rvindex)%dims(1) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'interpolate_nlat'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 1
    restartvars(rvindex)%dims(1) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'interpolate_nlon'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 1
    restartvars(rvindex)%dims(1) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'meridional_complement'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 2
    restartvars(rvindex)%dims(1) = maxnflds_dim_ind
    restartvars(rvindex)%dims(2) = ptapes_dim_ind

    rvindex = rvindex + 1
    restartvars(rvindex)%name = 'zonal_complement'
    restartvars(rvindex)%type = pio_int
    restartvars(rvindex)%ndims = 2
    restartvars(rvindex)%dims(1) = maxnflds_dim_ind
    restartvars(rvindex)%dims(2) = ptapes_dim_ind

  end subroutine restart_vars_setnames

  subroutine restart_dims_setnames()
    use cam_grid_support,    only: max_hcoordname_len
    use cam_history_support, only: registeredmdims

    restartdims(ptapes_dim_ind)%name = 'ptapes'
    restartdims(ptapes_dim_ind)%len  = ptapes

    restartdims(max_string_len_dim_ind)%name = 'max_string_len'
    restartdims(max_string_len_dim_ind)%len  = max_string_len

    restartdims(fieldname_lenp2_dim_ind)%name = 'fieldname_lenp2'
    restartdims(fieldname_lenp2_dim_ind)%len  = fieldname_lenp2

    restartdims(pflds_dim_ind)%name = 'pflds'
    restartdims(pflds_dim_ind)%len  = pflds

    restartdims(max_chars_dim_ind)%name = 'max_chars'
    restartdims(max_chars_dim_ind)%len  = max_chars

    restartdims(max_fieldname_len_dim_ind)%name = 'max_fieldname_len'
    restartdims(max_fieldname_len_dim_ind)%len  = max_fieldname_len

    restartdims(maxnflds_dim_ind)%name = 'maxnflds'
    restartdims(maxnflds_dim_ind)%len  = maxval(nflds)

    restartdims(maxvarmdims_dim_ind)%name = 'maxvarmdims'
    restartdims(maxvarmdims_dim_ind)%len  = maxvarmdims

    restartdims(registeredmdims_dim_ind)%name = 'registeredmdims'
    restartdims(registeredmdims_dim_ind)%len  = registeredmdims

    restartdims(max_hcoordname_len_dim_ind)%name = 'max_hcoordname_len'
    restartdims(max_hcoordname_len_dim_ind)%len  = max_hcoordname_len

  end subroutine restart_dims_setnames


  subroutine init_restart_history (File)
    use cam_pio_utils,  only: cam_pio_def_dim
    use cam_pio_utils,  only: cam_pio_handle_error

    !---------------------------------------------------------------------------
    !
    ! Arguments
    !
    type(file_desc_t), intent(inout) :: File                 ! Pio file Handle
    !
    ! Local
    !
    integer :: dimids(4), ndims
    integer :: ierr, i, k

    ! Don't need to write restart data if we have written the file this step
    where (hstwr(:))
      rgnht(:) = .false.
    elsewhere
      rgnht(:) = .true.
    end where

    if(maxval(nflds)>0) then
      call restart_vars_setnames()
      call restart_dims_setnames()

      do i=1,restartdimcnt
        ! it's possible that one or more of these have been defined elsewhere
        call cam_pio_def_dim(File, restartdims(i)%name, restartdims(i)%len,   &
             restartdims(i)%dimid, existOK=.true.)
      end do

      do i=1,restartvarcnt
        ndims= restartvars(i)%ndims
        do k=1,ndims
          dimids(k)=restartdims(restartvars(i)%dims(k))%dimid
        end do
        allocate(restartvars(i)%vdesc)
        ierr = pio_def_var(File, restartvars(i)%name, restartvars(i)%type, dimids(1:ndims), restartvars(i)%vdesc)
        call cam_pio_handle_error(ierr, 'INIT_RESTART_HISTORY: Error defining '//trim(restartvars(i)%name))

      end do
    end if
  end subroutine init_restart_history

  function restartvar_getdesc(name) result(vdesc)
    character(len=*), intent(in) :: name
    type(var_desc_t), pointer :: vdesc
    character(len=max_chars) :: errmsg
    integer :: i

    nullify(vdesc)
    do i=1,restartvarcnt
      if(name .eq. restartvars(i)%name) then
        vdesc=>restartvars(i)%vdesc
        exit
      end if
    end do
    if(.not.associated(vdesc)) then
      errmsg = 'Could not find restart variable '//name
      call endrun(errmsg)
    end if
  end function restartvar_getdesc


  !#######################################################################

  subroutine write_restart_history ( File, &
       yr_spec, mon_spec, day_spec, sec_spec )
    use cam_history_support, only: hist_coord_name, registeredmdims

    implicit none
    !--------------------------------------------------------------------------------------------------
    !
    ! Arguments
    !
    type(file_desc_t), intent(inout) :: file         ! PIO restart file pointer
    integer, intent(in), optional :: yr_spec         ! Simulation year
    integer, intent(in), optional :: mon_spec        ! Simulation month
    integer, intent(in), optional :: day_spec        ! Simulation day
    integer, intent(in), optional :: sec_spec        ! Seconds into current simulation day
    !
    ! Local workspace
    !
    integer :: ierr, t, f
    integer :: rgnht_int(ptapes), start(2), startc(3)
    type(var_desc_t), pointer :: vdesc

    ! PIO variable descriptors
    type(var_desc_t), pointer ::  field_name_desc   ! Restart field names
    type(var_desc_t), pointer ::  decomp_type_desc
    type(var_desc_t), pointer ::  numlev_desc
    type(var_desc_t), pointer ::  avgflag_desc
    type(var_desc_t), pointer ::  sseq_desc
    type(var_desc_t), pointer ::  cm_desc
    type(var_desc_t), pointer ::  longname_desc
    type(var_desc_t), pointer ::  units_desc
    type(var_desc_t), pointer ::  hwrt_prec_desc
    type(var_desc_t), pointer ::  xyfill_desc
    type(var_desc_t), pointer ::  mdims_desc        ! mdim name indices
    type(var_desc_t), pointer ::  mdimname_desc     ! mdim names
    type(var_desc_t), pointer ::  issubcol_desc
    type(var_desc_t), pointer ::  fillval_desc
    type(var_desc_t), pointer ::  interpolate_output_desc
    type(var_desc_t), pointer ::  interpolate_type_desc
    type(var_desc_t), pointer ::  interpolate_gridtype_desc
    type(var_desc_t), pointer ::  interpolate_nlat_desc
    type(var_desc_t), pointer ::  interpolate_nlon_desc
    type(var_desc_t), pointer ::  meridional_complement_desc
    type(var_desc_t), pointer ::  zonal_complement_desc

    integer, allocatable      ::  allmdims(:,:,:)
    integer, allocatable      ::  xyfill(:,:)
    integer, allocatable      ::  is_subcol(:,:)
    integer, allocatable      ::  interp_output(:)

    integer                   ::  maxnflds


    maxnflds = maxval(nflds)
    allocate(xyfill(maxnflds, ptapes))
    xyfill = 0
    allocate(is_subcol(maxnflds, ptapes))
    is_subcol = 0
    allocate(interp_output(ptapes))
    interp_output = 0

    !
    !-----------------------------------------------------------------------
    ! Write the history restart data if necessary
    !-----------------------------------------------------------------------

    rgnht_int(:) = 0

    if(.not.allocated(restarthistory_tape)) allocate(restarthistory_tape(ptapes))

    do t=1,ptapes
      ! No need to write history IC restart because it is always instantaneous
      if (is_initfile(file_index=t)) rgnht(t) = .false.
      ! No need to write restart data for empty files
      if (nflds(t) == 0) rgnht(t) = .false.
      if(rgnht(t)) then
        rgnht_int(t) = 1
        restarthistory_tape(t)%hlist => history_tape(t)%hlist

        if(associated(history_tape(t)%grid_ids)) then
          restarthistory_tape(t)%grid_ids => history_tape(t)%grid_ids
        end if
        if(associated(history_tape(t)%patches)) then
          restarthistory_tape(t)%patches => history_tape(t)%patches
        end if
      end if
    end do

    if(maxval(nflds)<=0) return

    call wshist(rgnht)

    vdesc => restartvar_getdesc('fincl')
    ierr= pio_put_var(File, vdesc, fincl(:,1:ptapes))

    vdesc => restartvar_getdesc('fincllonlat')
    ierr= pio_put_var(File, vdesc, fincllonlat(:,1:ptapes))

    vdesc => restartvar_getdesc('fexcl')
    ierr= pio_put_var(File, vdesc, fexcl(:,1:ptapes))

    vdesc => restartvar_getdesc('rgnht')
    ierr= pio_put_var(File, vdesc, rgnht_int(1:ptapes))

    vdesc => restartvar_getdesc('nhtfrq')
    ierr= pio_put_var(File, vdesc, nhtfrq(1:ptapes))

    vdesc => restartvar_getdesc('nflds')
    ierr= pio_put_var(File, vdesc, nflds(1:ptapes))

    vdesc => restartvar_getdesc('nfils')
    ierr= pio_put_var(File, vdesc, nfils(1:ptapes))

    vdesc => restartvar_getdesc('mfilt')
    ierr= pio_put_var(File, vdesc, mfilt(1:ptapes))

    vdesc => restartvar_getdesc('nfpath')
    ierr= pio_put_var(File, vdesc, nfpath(1:ptapes))

    vdesc => restartvar_getdesc('cpath')
    ierr= pio_put_var(File, vdesc,  cpath(1:ptapes))

    vdesc => restartvar_getdesc('nhfil')
    ierr= pio_put_var(File, vdesc, nhfil(1:ptapes))

    vdesc => restartvar_getdesc('ndens')
    ierr= pio_put_var(File, vdesc, ndens(1:ptapes))
    vdesc => restartvar_getdesc('ncprec')
    ierr= pio_put_var(File, vdesc, ncprec(1:ptapes))
    vdesc => restartvar_getdesc('beg_time')
    ierr= pio_put_var(File, vdesc, beg_time(1:ptapes))

    vdesc => restartvar_getdesc('hrestpath')
    ierr = pio_put_var(File, vdesc, hrestpath(1:ptapes))

    vdesc => restartvar_getdesc('lcltod_start')
    ierr = pio_put_var(File, vdesc, lcltod_start(1:ptapes))

    vdesc => restartvar_getdesc('lcltod_stop')
    ierr = pio_put_var(File, vdesc, lcltod_stop(1:ptapes))

    field_name_desc => restartvar_getdesc('field_name')
    decomp_type_desc => restartvar_getdesc('decomp_type')
    numlev_desc => restartvar_getdesc('numlev')
    hwrt_prec_desc => restartvar_getdesc('hwrt_prec')

    sseq_desc => restartvar_getdesc('sampling_seq')
    cm_desc => restartvar_getdesc('cell_methods')
    longname_desc => restartvar_getdesc('long_name')
    units_desc => restartvar_getdesc('units')
    avgflag_desc => restartvar_getdesc('avgflag')
    xyfill_desc => restartvar_getdesc('xyfill')
    issubcol_desc => restartvar_getdesc('is_subcol')

    interpolate_output_desc => restartvar_getdesc('interpolate_output')
    interpolate_type_desc => restartvar_getdesc('interpolate_type')
    interpolate_gridtype_desc => restartvar_getdesc('interpolate_gridtype')
    interpolate_nlat_desc => restartvar_getdesc('interpolate_nlat')
    interpolate_nlon_desc => restartvar_getdesc('interpolate_nlon')

    meridional_complement_desc => restartvar_getdesc('meridional_complement')
    zonal_complement_desc => restartvar_getdesc('zonal_complement')

    mdims_desc => restartvar_getdesc('mdims')
    mdimname_desc => restartvar_getdesc('mdimnames')
    fillval_desc => restartvar_getdesc('fillvalue')

    tape=>history_tape

    ! allmdims specifies the mdim indices for each field
    allocate(allmdims(maxvarmdims,maxval(nflds),ptapes))
    allmdims=-1

    startc(1)=1
    do t = 1,ptapes
      start(2)=t
      startc(3)=t
      do f=1,nflds(t)
        start(1)=f
        startc(2)=f
        ierr = pio_put_var(File, field_name_desc,startc,tape(t)%hlist(f)%field%name)
        ierr = pio_put_var(File, decomp_type_desc,start,tape(t)%hlist(f)%field%decomp_type)
        ierr = pio_put_var(File, numlev_desc,start,tape(t)%hlist(f)%field%numlev)

        ierr = pio_put_var(File, hwrt_prec_desc,start,tape(t)%hlist(f)%hwrt_prec)
        ierr = pio_put_var(File, sseq_desc,startc,tape(t)%hlist(f)%field%sampling_seq)
        ierr = pio_put_var(File, cm_desc,startc,tape(t)%hlist(f)%field%cell_methods)
        ierr = pio_put_var(File, longname_desc,startc,tape(t)%hlist(f)%field%long_name)
        ierr = pio_put_var(File, units_desc,startc,tape(t)%hlist(f)%field%units)
        ierr = pio_put_var(File, avgflag_desc,start, tape(t)%hlist(f)%avgflag)

        ierr = pio_put_var(File, fillval_desc,start, tape(t)%hlist(f)%field%fillvalue)
        ierr = pio_put_var(File, meridional_complement_desc,start, tape(t)%hlist(f)%field%meridional_complement)
        ierr = pio_put_var(File, zonal_complement_desc,start, tape(t)%hlist(f)%field%zonal_complement)
        if(associated(tape(t)%hlist(f)%field%mdims)) then
          allmdims(1:size(tape(t)%hlist(f)%field%mdims),f,t) = tape(t)%hlist(f)%field%mdims
        else
        end if
        if(tape(t)%hlist(f)%field%flag_xyfill) then
           xyfill(f,t) = 1
        end if
        if(tape(t)%hlist(f)%field%is_subcol) then
           is_subcol(f,t) = 1
        end if
      end do
      if (interpolate_output(t)) then
        interp_output(t) = 1
      end if
    end do
    ierr = pio_put_var(File, xyfill_desc, xyfill)
    ierr = pio_put_var(File, mdims_desc, allmdims)
    ierr = pio_put_var(File, issubcol_desc, is_subcol)
    !! Interpolated output variables
    ierr = pio_put_var(File, interpolate_output_desc, interp_output)
    interp_output = 1
    do t = 1, size(interpolate_info)
      interp_output(t) = interpolate_info(t)%interp_type
    end do
    ierr = pio_put_var(File, interpolate_type_desc, interp_output)
    interp_output = 1
    do t = 1, size(interpolate_info)
      interp_output(t) = interpolate_info(t)%interp_gridtype
    end do
    ierr = pio_put_var(File, interpolate_gridtype_desc, interp_output)
    interp_output = 0
    do t = 1, size(interpolate_info)
      interp_output(t) = interpolate_info(t)%interp_nlat
    end do
    ierr = pio_put_var(File, interpolate_nlat_desc, interp_output)
    interp_output = 0
    do t = 1, size(interpolate_info)
      interp_output(t) = interpolate_info(t)%interp_nlon
    end do
    ierr = pio_put_var(File, interpolate_nlon_desc, interp_output)
    ! Registered history coordinates
    start(1) = 1
    do f = 1, registeredmdims
      start(2) = f
      ierr = pio_put_var(File, mdimname_desc, start, hist_coord_name(f))
    end do

    deallocate(xyfill, allmdims)
    return

  end subroutine write_restart_history


  !#######################################################################

  subroutine read_restart_history (File)
    use pio,                 only: pio_inq_dimid
    use pio,                 only: pio_inq_varid, pio_inq_dimname
    use cam_pio_utils,       only: cam_pio_openfile, cam_pio_closefile
    use cam_pio_utils,       only: cam_pio_var_info
    use ioFileMod,           only: getfil
    use sat_hist,            only: sat_hist_define, sat_hist_init
    use cam_grid_support,    only: cam_grid_read_dist_array, cam_grid_num_grids
    use cam_history_support, only: get_hist_coord_index, add_hist_coord
    use constituents,        only: cnst_get_ind, cnst_get_type_byind

    use shr_sys_mod,         only: shr_sys_getenv
    use spmd_utils,          only: mpicom, mpi_character, masterprocid
    !
    !-----------------------------------------------------------------------
    !
    ! Arguments
    !
    type(file_desc_t), intent(inout) :: File            ! unit number
    !
    ! Local workspace
    !
    integer t, f, ff                 ! tape, field indices
    integer begdim2                  ! on-node vert start index
    integer enddim2                  ! on-node vert end index
    integer begdim1                  ! on-node dim1 start index
    integer enddim1                  ! on-node dim1 end index
    integer begdim3                  ! on-node chunk or lat start index
    integer enddim3                  ! on-node chunk or lat end index


    integer rgnht_int(ptapes)
    integer :: ierr

    character(len=max_string_len)  :: locfn       ! Local filename
    character(len=max_fieldname_len), allocatable :: tmpname(:,:)
    integer, allocatable :: decomp(:,:), tmpnumlev(:,:)
    integer, pointer :: nacs(:,:)    ! accumulation counter
    character(len=max_fieldname_len) :: fname_tmp ! local copy of field name
    character(len=max_fieldname_len) :: dname_tmp ! local copy of dim name

    integer :: i, ptapes_dimid

    type(var_desc_t)                 :: vdesc
    type(var_desc_t)                 :: longname_desc
    type(var_desc_t)                 :: units_desc
    type(var_desc_t)                 :: avgflag_desc
    type(var_desc_t)                 :: sseq_desc
    type(var_desc_t)                 :: cm_desc
    type(var_desc_t)                 :: fillval_desc
    type(var_desc_t)                 :: meridional_complement_desc
    type(var_desc_t)                 :: zonal_complement_desc
    integer,            allocatable  :: tmpprec(:,:)
    integer,            allocatable  :: xyfill(:,:)
    integer,            allocatable  :: allmdims(:,:,:)
    integer,            allocatable  :: is_subcol(:,:)
    integer,            allocatable  :: interp_output(:)
    integer                          :: nacsdimcnt, nacsval
    integer                          :: maxnflds, dimid

    ! List of active grids (first dim) for each tape (second dim)
    ! An active grid is one for which there is a least one field being output
    !    on that grid.
    integer, allocatable        :: gridsontape(:,:)

    character(len=16),  allocatable  :: mdimnames(:) ! Names of all hist coords (inc. vertical)
    integer                          :: ndims, dimids(8)
    integer                          :: tmpdims(8), dimcnt
    integer                          :: dimlens(7)
    integer                          :: mtapes, mdimcnt
    integer                          :: fdims(3)         ! Field dims
    integer                          :: nfdims           ! 2 or 3 (for 2D,3D)
    integer                          :: fdecomp          ! Grid ID for field
    integer                          :: idx
    character(len=3)                 :: mixing_ratio

    !
    ! Get users logname and machine hostname
    !
    if ( masterproc )then
      logname = ' '
      call shr_sys_getenv ('LOGNAME',logname,ierr)
      host = ' '
      call shr_sys_getenv ('HOST',host,ierr)
    end if
    ! PIO requires netcdf attributes have consistant values on all tasks
    call mpi_bcast(logname, len(logname), mpi_character, masterprocid, mpicom, ierr)
    call mpi_bcast(host,    len(host),    mpi_character, masterprocid, mpicom, ierr)

    call pio_seterrorhandling(File, PIO_BCAST_ERROR)

    ierr = pio_inq_dimid(File, 'ptapes', ptapes_dimid)
    if(ierr/= PIO_NOERR) then
      if(masterproc) write(iulog,*) 'Not reading history info from restart file', ierr
      return   ! no history info in restart file
    end if
    call pio_seterrorhandling(File, PIO_INTERNAL_ERROR)

    ierr = pio_inq_dimlen(File, ptapes_dimid, mtapes)

    ierr = pio_inq_dimid(File, 'maxnflds', dimid)
    ierr = pio_inq_dimlen(File, dimid, maxnflds)

    ierr = pio_inq_dimid(File, 'maxvarmdims', dimid)
    ierr = pio_inq_dimlen(File, dimid, maxvarmdims)

    ierr = pio_inq_varid(File, 'rgnht', vdesc)
    ierr = pio_get_var(File, vdesc, rgnht_int(1:mtapes))

    ierr = pio_inq_varid(File, 'nhtfrq', vdesc)
    ierr = pio_get_var(File, vdesc, nhtfrq(1:mtapes))

    ierr = pio_inq_varid(File, 'nflds', vdesc)
    ierr = pio_get_var(File, vdesc, nflds(1:mtapes))
    ierr = pio_inq_varid(File, 'nfils', vdesc)
    ierr = pio_get_var(File, vdesc, nfils(1:mtapes))
    ierr = pio_inq_varid(File, 'mfilt', vdesc)
    ierr = pio_get_var(File, vdesc, mfilt(1:mtapes))

    ierr = pio_inq_varid(File, 'nfpath', vdesc)
    ierr = pio_get_var(File, vdesc, nfpath(1:mtapes))
    ierr = pio_inq_varid(File, 'cpath', vdesc)
    ierr = pio_get_var(File, vdesc, cpath(1:mtapes))
    ierr = pio_inq_varid(File, 'nhfil', vdesc)
    ierr = pio_get_var(File, vdesc, nhfil(1:mtapes))
    ierr = pio_inq_varid(File, 'hrestpath', vdesc)
    ierr = pio_get_var(File, vdesc, hrestpath(1:mtapes))


    ierr = pio_inq_varid(File, 'ndens', vdesc)
    ierr = pio_get_var(File, vdesc, ndens(1:mtapes))
    ierr = pio_inq_varid(File, 'ncprec', vdesc)
    ierr = pio_get_var(File, vdesc, ncprec(1:mtapes))
    ierr = pio_inq_varid(File, 'beg_time', vdesc)
    ierr = pio_get_var(File, vdesc, beg_time(1:mtapes))


    ierr = pio_inq_varid(File, 'fincl', vdesc)
    ierr = pio_get_var(File, vdesc, fincl(:,1:mtapes))

    ierr = pio_inq_varid(File, 'fincllonlat', vdesc)
    ierr = pio_get_var(File, vdesc, fincllonlat(:,1:mtapes))

    ierr = pio_inq_varid(File, 'fexcl', vdesc)
    ierr = pio_get_var(File, vdesc, fexcl(:,1:mtapes))

    ierr = pio_inq_varid(File, 'lcltod_start', vdesc)
    ierr = pio_get_var(File, vdesc, lcltod_start(1:mtapes))

    ierr = pio_inq_varid(File, 'lcltod_stop', vdesc)
    ierr = pio_get_var(File, vdesc, lcltod_stop(1:mtapes))




    allocate(tmpname(maxnflds, mtapes), decomp(maxnflds, mtapes), tmpnumlev(maxnflds,mtapes))
    ierr = pio_inq_varid(File, 'field_name', vdesc)
    ierr = pio_get_var(File, vdesc, tmpname)

    ierr = pio_inq_varid(File, 'decomp_type', vdesc)
    ierr = pio_get_var(File, vdesc, decomp)
    ierr = pio_inq_varid(File, 'numlev', vdesc)
    ierr = pio_get_var(File, vdesc, tmpnumlev)

    allocate(tmpprec(maxnflds,mtapes))
    ierr = pio_inq_varid(File, 'hwrt_prec',vdesc)
    ierr = pio_get_var(File, vdesc, tmpprec(:,:))

    allocate(xyfill(maxnflds,mtapes))
    ierr = pio_inq_varid(File, 'xyfill', vdesc)
    ierr = pio_get_var(File, vdesc, xyfill)

    allocate(is_subcol(maxnflds,mtapes))
    ierr = pio_inq_varid(File, 'is_subcol', vdesc)
    ierr = pio_get_var(File, vdesc, is_subcol)

    !! interpolated output
    allocate(interp_output(mtapes))
    ierr = pio_inq_varid(File, 'interpolate_output', vdesc)
    ierr = pio_get_var(File, vdesc, interp_output)
    interpolate_output(1:mtapes) = interp_output(1:mtapes) > 0
    if (ptapes > mtapes) then
      interpolate_output(mtapes+1:ptapes) = .false.
    end if
    ierr = pio_inq_varid(File, 'interpolate_type', vdesc)
    ierr = pio_get_var(File, vdesc, interp_output)
    do t = 1, mtapes
      if (interpolate_output(t)) then
        interpolate_info(t)%interp_type = interp_output(t)
      end if
    end do
    ierr = pio_inq_varid(File, 'interpolate_gridtype', vdesc)
    ierr = pio_get_var(File, vdesc, interp_output)
    do t = 1, mtapes
      if (interpolate_output(t)) then
        interpolate_info(t)%interp_gridtype = interp_output(t)
      end if
    end do
    ierr = pio_inq_varid(File, 'interpolate_nlat', vdesc)
    ierr = pio_get_var(File, vdesc, interp_output)
    do t = 1, mtapes
      if (interpolate_output(t)) then
        interpolate_info(t)%interp_nlat = interp_output(t)
      end if
    end do
    ierr = pio_inq_varid(File, 'interpolate_nlon', vdesc)
    ierr = pio_get_var(File, vdesc, interp_output)
    do t = 1, mtapes
      if (interpolate_output(t)) then
        interpolate_info(t)%interp_nlon = interp_output(t)
      end if
    end do

    !! mdim indices
    allocate(allmdims(maxvarmdims,maxnflds,mtapes))
    ierr = pio_inq_varid(File, 'mdims', vdesc)
    ierr = pio_get_var(File, vdesc, allmdims)

    !! mdim names
    ! Read the hist coord names to make sure they are all registered
    ierr = pio_inq_varid(File, 'mdimnames', vdesc)
    call cam_pio_var_info(File, vdesc, ndims, dimids, dimlens)
    mdimcnt = dimlens(2)
    allocate(mdimnames(mdimcnt))
    ierr = pio_get_var(File, vdesc, mdimnames)
    do f = 1, mdimcnt
      ! Check to see if the mdim is registered
      if (get_hist_coord_index(trim(mdimnames(f))) <= 0) then
        ! We need to register this mdim (hist_coord)
        call add_hist_coord(trim(mdimnames(f)))
      end if
    end do

    ierr = pio_inq_varid(File, 'avgflag', avgflag_desc)

    ierr = pio_inq_varid(File, 'long_name', longname_desc)
    ierr = pio_inq_varid(File, 'units', units_desc)
    ierr = pio_inq_varid(File, 'sampling_seq', sseq_desc)
    ierr = pio_inq_varid(File, 'cell_methods', cm_desc)

    ierr = pio_inq_varid(File, 'fillvalue', fillval_desc)
    ierr = pio_inq_varid(File, 'meridional_complement', meridional_complement_desc)
    ierr = pio_inq_varid(File, 'zonal_complement', zonal_complement_desc)

    rgnht(:)=.false.

    allocate(history_tape(mtapes))

    tape => history_tape

    do t=1,mtapes

      if(rgnht_int(t)==1) rgnht(t)=.true.


      call strip_null(nfpath(t))
      call strip_null(cpath(t))
      call strip_null(hrestpath(t))
      allocate(tape(t)%hlist(nflds(t)))

      do f=1,nflds(t)
        if (associated(tape(t)%hlist(f)%field%mdims)) then
          deallocate(tape(t)%hlist(f)%field%mdims)
        end if
        nullify(tape(t)%hlist(f)%field%mdims)
        ierr = pio_get_var(File,fillval_desc, (/f,t/), tape(t)%hlist(f)%field%fillvalue)
        ierr = pio_get_var(File,meridional_complement_desc, (/f,t/), tape(t)%hlist(f)%field%meridional_complement)
        ierr = pio_get_var(File,zonal_complement_desc, (/f,t/), tape(t)%hlist(f)%field%zonal_complement)
        ierr = pio_get_var(File,avgflag_desc, (/f,t/), tape(t)%hlist(f)%avgflag)
        ierr = pio_get_var(File,longname_desc, (/1,f,t/), tape(t)%hlist(f)%field%long_name)
        ierr = pio_get_var(File,units_desc, (/1,f,t/), tape(t)%hlist(f)%field%units)
        tape(t)%hlist(f)%field%sampling_seq(1:max_chars) = ' '
        ierr = pio_get_var(File,sseq_desc, (/1,f,t/), tape(t)%hlist(f)%field%sampling_seq)
        call strip_null(tape(t)%hlist(f)%field%sampling_seq)
        tape(t)%hlist(f)%field%cell_methods(1:max_chars) = ' '
        ierr = pio_get_var(File,cm_desc, (/1,f,t/), tape(t)%hlist(f)%field%cell_methods)
        call strip_null(tape(t)%hlist(f)%field%cell_methods)
        if(xyfill(f,t) ==1) then
          tape(t)%hlist(f)%field%flag_xyfill=.true.
        else
          tape(t)%hlist(f)%field%flag_xyfill=.false.
        end if
        if(is_subcol(f,t) ==1) then
           tape(t)%hlist(f)%field%is_subcol=.true.
        else
           tape(t)%hlist(f)%field%is_subcol=.false.
        end if
        call strip_null(tmpname(f,t))
        tape(t)%hlist(f)%field%name = tmpname(f,t)
        tape(t)%hlist(f)%field%decomp_type = decomp(f,t)
        tape(t)%hlist(f)%field%numlev = tmpnumlev(f,t)
        tape(t)%hlist(f)%hwrt_prec = tmpprec(f,t)

        ! If the field is an advected constituent set the mixing_ratio attribute
        fname_tmp = strip_suffix(tape(t)%hlist(f)%field%name)
        call cnst_get_ind(fname_tmp, idx, abort=.false.)
        mixing_ratio = ''
        if (idx > 0) then
           mixing_ratio = cnst_get_type_byind(idx)
        end if
        tape(t)%hlist(f)%field%mixing_ratio = mixing_ratio

        mdimcnt = count(allmdims(:,f,t) > 0)
        if(mdimcnt > 0) then
          allocate(tape(t)%hlist(f)%field%mdims(mdimcnt))
          do i = 1, mdimcnt
            tape(t)%hlist(f)%field%mdims(i) = get_hist_coord_index(mdimnames(allmdims(i,f,t)))
          end do
        end if

      end do
    end do
    deallocate(tmpname, tmpnumlev, tmpprec, decomp, xyfill, is_subcol)
    deallocate(mdimnames)

    allocate(gridsontape(cam_grid_num_grids() + 1, ptapes))
    gridsontape = -1
    do t = 1, ptapes
      do f = 1, nflds(t)
        call set_field_dimensions(tape(t)%hlist(f)%field)

        begdim1 = tape(t)%hlist(f)%field%begdim1
        enddim1 = tape(t)%hlist(f)%field%enddim1
        begdim2 = tape(t)%hlist(f)%field%begdim2
        enddim2 = tape(t)%hlist(f)%field%enddim2
        begdim3 = tape(t)%hlist(f)%field%begdim3
        enddim3 = tape(t)%hlist(f)%field%enddim3

        allocate(tape(t)%hlist(f)%hbuf(begdim1:enddim1,begdim2:enddim2,begdim3:enddim3))
        if (tape(t)%hlist(f)%avgflag .eq. 'S') then ! allocate the variance buffer for standard dev
           allocate(tape(t)%hlist(f)%sbuf(begdim1:enddim1,begdim2:enddim2,begdim3:enddim3))
        endif

        if (associated(tape(t)%hlist(f)%varid)) then
          deallocate(tape(t)%hlist(f)%varid)
        end if
        nullify(tape(t)%hlist(f)%varid)
        if (associated(tape(t)%hlist(f)%nacs)) then
          deallocate(tape(t)%hlist(f)%nacs)
        end if
        nullify(tape(t)%hlist(f)%nacs)
        if(tape(t)%hlist(f)%field%flag_xyfill .or. (avgflag_pertape(t)=='L')) then
          allocate (tape(t)%hlist(f)%nacs(begdim1:enddim1,begdim3:enddim3))
        else
          allocate(tape(t)%hlist(f)%nacs(1,begdim3:enddim3))
        end if
        ! initialize all buffers to zero - this will be overwritten later by the
        ! data in the history restart file if it exists.
        call h_zero(f,t)

        ! Make sure this field's decomp is listed on the tape
        fdecomp = tape(t)%hlist(f)%field%decomp_type
        do ff = 1, size(gridsontape, 1)
          if (fdecomp == gridsontape(ff, t)) then
            exit
          else if (gridsontape(ff, t) < 0) then
            gridsontape(ff, t) = fdecomp
            exit
          end if
        end do

      end do
    end do
    !
    !-----------------------------------------------------------------------
    ! Read history restart files
    !-----------------------------------------------------------------------
    !
    ! Loop over the total number of history files declared and
    ! read the pathname for any history restart files
    ! that are present (if any). Test to see if the run is a restart run
    ! AND if any history buffer regen files exist (rgnht=.T.). Note, rgnht
    ! is preset to false, reset to true in routine WSDS if hbuf restart files
    ! are written and saved in the master restart file. Each history buffer
    ! restart file is then obtained.
    ! Note: some f90 compilers (e.g. SGI) complain about I/O of
    ! derived types which have pointer components, so explicitly read each one.
    !
    do t=1,mtapes
      if (rgnht(t)) then
        !
        ! Open history restart file
        !
        call getfil (hrestpath(t), locfn)
        call cam_pio_openfile(tape(t)%File, locfn, 0)
        !
        ! Read history restart file
        !
        do f = 1, nflds(t)

          fname_tmp = strip_suffix(tape(t)%hlist(f)%field%name)
          if(masterproc) write(iulog, *) 'Reading history variable ',fname_tmp
          ierr = pio_inq_varid(tape(t)%File, fname_tmp, vdesc)

          call cam_pio_var_info(tape(t)%File, vdesc, ndims, dimids, dimlens)
          if(.not. associated(tape(t)%hlist(f)%field%mdims)) then
            dimcnt = 0
            do i=1,ndims
              ierr = pio_inq_dimname(tape(t)%File, dimids(i), dname_tmp)
              dimid = get_hist_coord_index(dname_tmp)
              if(dimid >= 1) then
                dimcnt = dimcnt + 1
                tmpdims(dimcnt) = dimid
              ! No else, just looking for mdims (grid dims won't be hist coords)
              end if
            end do
            if(dimcnt > 0) then
              allocate(tape(t)%hlist(f)%field%mdims(dimcnt))
              tape(t)%hlist(f)%field%mdims(:) = tmpdims(1:dimcnt)
              if(dimcnt > maxvarmdims) maxvarmdims=dimcnt
            end if
          end if
          call set_field_dimensions(tape(t)%hlist(f)%field)
          begdim1    =  tape(t)%hlist(f)%field%begdim1
          enddim1    =  tape(t)%hlist(f)%field%enddim1
          fdims(1)   =  enddim1 - begdim1 + 1
          begdim2    =  tape(t)%hlist(f)%field%begdim2
          enddim2    =  tape(t)%hlist(f)%field%enddim2
          fdims(2)   =  enddim2 - begdim2 + 1
          begdim3    =  tape(t)%hlist(f)%field%begdim3
          enddim3    =  tape(t)%hlist(f)%field%enddim3
          fdims(3)   =  enddim3 - begdim3 + 1
          if (fdims(2) > 1) then
            nfdims = 3
          else
            nfdims = 2
            fdims(2) = fdims(3)
          end if
          fdecomp = tape(t)%hlist(f)%field%decomp_type
          if (nfdims > 2) then
            call cam_grid_read_dist_array(tape(t)%File, fdecomp,              &
                 fdims(1:nfdims), dimlens(1:ndims), tape(t)%hlist(f)%hbuf, vdesc)
          else
            call cam_grid_read_dist_array(tape(t)%File, fdecomp,              &
                 fdims(1:nfdims), dimlens(1:ndims), tape(t)%hlist(f)%hbuf(:,1,:), vdesc)
          end if

          if ( associated(tape(t)%hlist(f)%sbuf) ) then
             ! read in variance for standard deviation
             ierr = pio_inq_varid(tape(t)%File, trim(fname_tmp)//'_var', vdesc)
             if (nfdims > 2) then
                call cam_grid_read_dist_array(tape(t)%File, fdecomp,              &
                     fdims(1:nfdims), dimlens(1:ndims), tape(t)%hlist(f)%sbuf, vdesc)
             else
                call cam_grid_read_dist_array(tape(t)%File, fdecomp,              &
                     fdims(1:nfdims), dimlens(1:ndims), tape(t)%hlist(f)%sbuf(:,1,:), vdesc)
             end if
          endif

          ierr = pio_inq_varid(tape(t)%File, trim(fname_tmp)//'_nacs', vdesc)
          call cam_pio_var_info(tape(t)%File, vdesc, nacsdimcnt, dimids, dimlens)

          if(nacsdimcnt > 0) then
            if (nfdims > 2) then
              ! nacs only has 2 dims (no levels)
              fdims(2) = fdims(3)
            end if
            allocate(tape(t)%hlist(f)%nacs(begdim1:enddim1,begdim3:enddim3))
            nacs       => tape(t)%hlist(f)%nacs(:,:)
            call cam_grid_read_dist_array(tape(t)%File, fdecomp, fdims(1:2),  &
                 dimlens(1:nacsdimcnt), nacs, vdesc)
          else
            allocate(tape(t)%hlist(f)%nacs(1,begdim3:enddim3))
            ierr = pio_get_var(tape(t)%File, vdesc, nacsval)
            tape(t)%hlist(f)%nacs(1,:)= nacsval
          end if

        end do
        !
        ! Done reading this history restart file
        !
        call cam_pio_closefile(tape(t)%File)

      end if  ! rgnht(t)

      ! (re)create the master list of grid IDs
      ff = 0
      do f = 1, size(gridsontape, 1)
        if (gridsontape(f, t) > 0) then
          ff = ff + 1
        end if
      end do
      allocate(tape(t)%grid_ids(ff))
      ff = 1
      do f = 1, size(gridsontape, 1)
        if (gridsontape(f, t) > 0) then
          tape(t)%grid_ids(ff) = gridsontape(f, t)
          ff = ff + 1
        end if
      end do
      call patch_init(t)
    end do     ! end of do mtapes loop

    !
    ! If the history files are partially complete (contain less than
    ! mfilt(t) time samples, then get the files and open them.)
    !
    ! NOTE:  No need to perform this operation for IC history files or empty files
    !

    do t=1,mtapes
      if (is_initfile(file_index=t)) then
        ! Initialize filename specifier for IC file
        hfilename_spec(t) = '%c.cam' // trim(inst_suffix) // '.i.%y-%m-%d-%s.nc'
        nfils(t) = 0
      else if (nflds(t) == 0) then
        nfils(t) = 0
      else
        if (nfils(t) > 0) then
          call getfil (cpath(t), locfn)
          call cam_pio_openfile(tape(t)%File, locfn, PIO_WRITE)
          call h_inquire (t)
          if(is_satfile(t)) then
            !  Initialize the sat following history subsystem
            call sat_hist_init()
            call sat_hist_define(tape(t)%File)
          end if
        end if
        !
        ! If the history file is full, close the current unit
        !
        if (nfils(t) >= mfilt(t)) then
          if (masterproc) then
            write(iulog,*)'READ_RESTART_HISTORY: nf_close(',t,')=',nhfil(t), mfilt(t)
          end if
          do f=1,nflds(t)
            deallocate(tape(t)%hlist(f)%varid)
            nullify(tape(t)%hlist(f)%varid)
          end do
          call cam_pio_closefile(tape(t)%File)
          nfils(t) = 0
        end if
      end if
    end do

    ! Setup vector pairs for unstructured grid interpolation
    call setup_interpolation_and_define_vector_complements()

    if(mtapes/=ptapes .and. masterproc) then
      write(iulog,*) ' WARNING: Restart file ptapes setting ',mtapes,' not equal to model setting ',ptapes
    end if

    return
  end subroutine read_restart_history

  !#######################################################################

  character(len=max_string_len) function get_hfilepath( tape )
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Return full filepath of history file for given tape number
    ! This allows public read access to the filenames without making
    ! the filenames public data.
    !
    !-----------------------------------------------------------------------
    !
    integer, intent(in) :: tape  ! Tape number

    get_hfilepath = cpath( tape )
  end function get_hfilepath

  !#######################################################################

  character(len=max_string_len) function get_hist_restart_filepath( tape )
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Return full filepath of restart file for given tape number
    ! This allows public read access to the filenames without making
    ! the filenames public data.
    !
    !-----------------------------------------------------------------------
    !
    integer, intent(in) :: tape  ! Tape number

    get_hist_restart_filepath = hrestpath( tape )
  end function get_hist_restart_filepath

  !#######################################################################

  integer function get_ptapes( )
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Return the number of tapes being used.
    ! This allows public read access to the number of tapes without making
    ! ptapes public data.
    !
    !-----------------------------------------------------------------------
    !
    get_ptapes = ptapes
  end function get_ptapes

  !#######################################################################

  recursive function get_entry_by_name(listentry, name) result(entry)
    type(master_entry),  pointer :: listentry
    character(len=*), intent(in) :: name ! variable name
    type(master_entry), pointer :: entry

    if(associated(listentry)) then
      if(listentry%field%name .eq. name) then
        entry => listentry
      else
        entry=>get_entry_by_name(listentry%next_entry, name)
      end if
    else
      nullify(entry)
    end if
  end function get_entry_by_name

  !#######################################################################

  subroutine AvgflagToString(avgflag, time_op)
    ! Dummy arguments
    character(len=1),           intent(in)  :: avgflag ! averaging flag
    character(len=max_chars),   intent(out) :: time_op ! time op (e.g. max)

    ! Local variable
    character(len=*), parameter             :: subname = 'AvgflagToString'

    select case (avgflag)
    case ('A')
      time_op(:) = 'mean'
    case ('B')
      time_op(:) = 'mean00z'
    case ('I')
      time_op(:) = ' '
    case ('X')
      time_op(:) = 'maximum'
    case ('M')
      time_op(:) = 'minimum'
    case('L')
      time_op(:) = LT_DESC
    case ('S')
      time_op(:) = 'standard_deviation'
    case default
      call endrun(subname//': unknown avgflag = '//avgflag)
    end select
  end subroutine AvgflagToString

  !#######################################################################

  subroutine fldlst ()

    use cam_grid_support, only: cam_grid_num_grids
    use spmd_utils,       only: mpicom
    use dycore,           only: dycore_is

    !-----------------------------------------------------------------------
    !
    ! Purpose: Define the contents of each history file based on namelist input for initial or branch
    ! run, and restart data if a restart run.
    !
    ! Method: Use arrays fincl and fexcl to modify default history tape contents.
    !         Then sort the result alphanumerically for later use by OUTFLD to
    !         allow an n log n search time.
    !
    !---------------------------Local variables-----------------------------
    !
    integer t, f                   ! tape, field indices
    integer ff                     ! index into include, exclude and fprec list
    integer :: i
    character(len=fieldname_len) :: name ! field name portion of fincl (i.e. no avgflag separator)
    character(len=max_fieldname_len) :: mastername ! name from masterlist field
    character(len=max_chars) :: errormsg ! error output field
    character(len=1) :: avgflag    ! averaging flag
    character(len=1) :: prec_wrt   ! history buffer write precision flag

    type (hentry) :: tmp           ! temporary used for swapping

    type(master_entry), pointer :: listentry
    logical                     :: fieldontape      ! .true. iff field on tape
    integer                     :: errors_found

    ! List of active grids (first dim) for each tape (second dim)
    ! An active grid is one for which there is a least one field being output
    !    on that grid.
    integer, allocatable        :: gridsontape(:,:)

    ! The following list of field names are only valid for the FV dycore.  They appear
    ! in fincl settings of WACCM use case files which are not restricted to the FV dycore.
    ! To avoid duplicating long fincl lists in use case files to provide both FV and non-FV
    ! versions this short list of fields is checked for and removed from fincl lists when
    ! the dycore is not FV.
    integer, parameter :: n_fv_only = 10
    character(len=6) :: fv_only_flds(n_fv_only) = &
       [ 'VTHzm ', 'WTHzm ', 'UVzm  ', 'UWzm  ', 'Uzm   ', 'Vzm   ', 'Wzm   ', &
         'THzm  ', 'TH    ', 'MSKtem' ]

    integer :: n_vec_comp, add_fincl_idx
    integer, parameter :: nvecmax = 50 ! max number of vector components in a fincl list
    character(len=2) :: avg_suffix
    character(len=max_fieldname_len) :: vec_comp_names(nvecmax)
    character(len=1)                 :: vec_comp_avgflag(nvecmax)
    !--------------------------------------------------------------------------

    ! First ensure contents of fincl, fexcl, and fwrtpr are all valid names
    !
    errors_found = 0
    do t=1,ptapes

      f = 1
      n_vec_comp       = 0
      vec_comp_names   = ' '
      vec_comp_avgflag = ' '
fincls: do while (f < pflds .and. fincl(f,t) /= ' ')
        name = getname (fincl(f,t))

        if (.not. dycore_is('FV')) then
           ! filter out fields only provided by FV dycore
           do i = 1, n_fv_only
              if (name == fv_only_flds(i)) then
                 write(errormsg,'(3a,2(i0,a))')'FLDLST: ', trim(name), &
                    ' in fincl(', f,', ',t, ') only available with FV dycore'
                 if (masterproc) then
                    write(iulog,*) trim(errormsg)
                    call shr_sys_flush(iulog)
                 end if
                 f = f + 1
                 cycle fincls
              end if
           end do
        end if

        mastername=''
        listentry => get_entry_by_name(masterlinkedlist, name)
        if (associated(listentry)) mastername = listentry%field%name
        if (name /= mastername) then
          write(errormsg,'(3a,2(i0,a))')'FLDLST: ', trim(name), ' in fincl(', f,', ',t, ') not found'
          if (masterproc) then
             write(iulog,*) trim(errormsg)
             call shr_sys_flush(iulog)
          end if
          errors_found = errors_found + 1
        else
           if (len_trim(mastername)>0 .and. interpolate_output(t)) then
              if (n_vec_comp >= nvecmax) call endrun('FLDLST: need to increase nvecmax')
              ! If this is a vector component then save the name of the complement
              avgflag = getflag(fincl(f,t))
              if (len_trim(listentry%meridional_field) > 0) then
                 n_vec_comp = n_vec_comp + 1
                 vec_comp_names(n_vec_comp) = listentry%meridional_field
                 vec_comp_avgflag(n_vec_comp) = avgflag
              else if (len_trim(listentry%zonal_field) > 0) then
                 n_vec_comp = n_vec_comp + 1
                 vec_comp_names(n_vec_comp) = listentry%zonal_field
                 vec_comp_avgflag(n_vec_comp) = avgflag
              end if
           end if
        end if
        f = f + 1
      end do fincls

      ! Interpolation of vector components requires that both be present.  If the fincl
      ! specifier contains any vector components, then the complement was saved in the
      ! array vec_comp_names.  Next insure (for interpolated output only) that all complements
      ! are also present in the fincl array.

      ! The first empty slot in the current fincl array is index f from loop above.
      add_fincl_idx = f
      if (f > 1 .and. interpolate_output(t)) then
         do i = 1, n_vec_comp
            call list_index(fincl(:,t), vec_comp_names(i), ff)
            if (ff == 0) then

               ! Add vector component to fincl.  Don't need to check whether its in the master
               ! list since this was done at the time of registering the vector components.
               avg_suffix = '  '
               if (len_trim(vec_comp_avgflag(i)) > 0) avg_suffix = ':' // vec_comp_avgflag(i)
               fincl(add_fincl_idx,t) = trim(vec_comp_names(i)) // avg_suffix
               add_fincl_idx = add_fincl_idx + 1

               write(errormsg,'(3a,1(i0,a))')'FLDLST: ', trim(vec_comp_names(i)), &
                  ' added to fincl', t, '.  Both vector components are required for interpolated output.'
               if (masterproc) then
                  write(iulog,*) trim(errormsg)
                  call shr_sys_flush(iulog)
               end if
            end if
         end do
      end if

      f = 1
      do while (f < pflds .and. fexcl(f,t) /= ' ')
        mastername=''
        listentry => get_entry_by_name(masterlinkedlist, fexcl(f,t))
        if(associated(listentry)) mastername = listentry%field%name

        if (fexcl(f,t) /= mastername) then
          write(errormsg,'(3a,2(i0,a))')'FLDLST: ', trim(fexcl(f,t)), ' in fexcl(', f,', ',t, ') not found'
          if (masterproc) then
             write(iulog,*) trim(errormsg)
             call shr_sys_flush(iulog)
          end if
          errors_found = errors_found + 1
        end if
        f = f + 1
      end do

      f = 1
      do while (f < pflds .and. fwrtpr(f,t) /= ' ')
        name = getname (fwrtpr(f,t))
        mastername=''
        listentry => get_entry_by_name(masterlinkedlist, name)
        if(associated(listentry)) mastername = listentry%field%name
        if (name /= mastername) then
          write(errormsg,'(3a,i0,a)')'FLDLST: ', trim(name), ' in fwrtpr(', f, ') not found'
          if (masterproc) then
             write(iulog,*) trim(errormsg)
             call shr_sys_flush(iulog)
          end if
          errors_found = errors_found + 1
        end if
        do ff=1,f-1                 ! If duplicate entry is found, stop
          if (trim(name) == trim(getname(fwrtpr(ff,t)))) then
            write(errormsg,'(3a)')'FLDLST: Duplicate field ', trim(name), ' in fwrtpr'
            if (masterproc) then
               write(iulog,*) trim(errormsg)
               call shr_sys_flush(iulog)
            end if
            errors_found = errors_found + 1
          end if
        end do
        f = f + 1
      end do
    end do

    if (errors_found > 0) then
       ! Give masterproc a chance to write all the log messages
       call mpi_barrier(mpicom, t)
       write(errormsg, '(a,i0,a)') 'FLDLST: ',errors_found,' errors found, see log'
       call endrun(trim(errormsg))
    end if

    nflds(:) = 0
    ! IC history file is to be created, set properties
    if(is_initfile()) then
      hfilename_spec(ptapes) = '%c.cam' // trim(inst_suffix) // '.i.%y-%m-%d-%s.nc'

      ncprec(ptapes) = pio_double
      ndens (ptapes) = 1
      mfilt (ptapes) = 1
    end if



    allocate(gridsontape(cam_grid_num_grids() + 1, ptapes))
    gridsontape = -1
    do t=1,ptapes
      !
      ! Add the field to the tape if specified via namelist (FINCL[1-ptapes]), or if
      ! it is on by default and was not excluded via namelist (FEXCL[1-ptapes]).
      ! Also set history buffer accumulation and output precision values according
      ! to the values specified via namelist (FWRTPR[1-ptapes])
      ! or, if not on the list, to the default values given by ndens(t).
      !
      listentry => masterlinkedlist
      do while(associated(listentry))
        mastername = listentry%field%name
        call list_index (fincl(1,t), mastername, ff)

        fieldontape = .false.
        if (ff > 0) then
          fieldontape = .true.
        else if ((.not. empty_htapes) .or. (is_initfile(file_index=t))) then
          call list_index (fexcl(1,t), mastername, ff)
          if (ff == 0 .and. listentry%actflag(t)) then
            fieldontape = .true.
          end if
        end if
        if (fieldontape) then
          ! The field is active so increment the number fo fields and add
          ! its decomp type to the list of decomp types on this tape
          nflds(t) = nflds(t) + 1
          do ff = 1, size(gridsontape, 1)
            if (listentry%field%decomp_type == gridsontape(ff, t)) then
              exit
            else if (gridsontape(ff, t) < 0) then
              gridsontape(ff, t) = listentry%field%decomp_type
              exit
            end if
          end do
        end if
        listentry=>listentry%next_entry
      end do
    end do
    !
    ! Determine total number of active history tapes
    !
    if (masterproc) then
      do t=1,ptapes
        if (nflds(t)  ==  0) then
          write(iulog,*)'FLDLST: Tape ',t,' is empty'
        end if
      end do
    endif
    allocate(history_tape(ptapes))
    tape=>history_tape


    do t=1,ptapes
      nullify(tape(t)%hlist)
      ! Now we have a field count and can allocate
      if(nflds(t) > 0) then
        ! Allocate the correct number of hentry slots
        allocate(tape(t)%hlist(nflds(t)))
        ! Count up the number of grids output on this tape
        ff = 0
        do f = 1, size(gridsontape, 1)
          if (gridsontape(f, t) > 0) then
            ff = ff + 1
          end if
        end do
        allocate(tape(t)%grid_ids(ff))
        ff = 1
        do f = 1, size(gridsontape, 1)
          if (gridsontape(f, t) > 0) then
            tape(t)%grid_ids(ff) = gridsontape(f, t)
            ff = ff + 1
          end if
        end do
      end if
      do ff=1,nflds(t)
        nullify(tape(t)%hlist(ff)%hbuf)
        nullify(tape(t)%hlist(ff)%sbuf)
        nullify(tape(t)%hlist(ff)%nacs)
        nullify(tape(t)%hlist(ff)%varid)
      end do


      nflds(t) = 0 ! recount to support array based method
      listentry => masterlinkedlist
      do while(associated(listentry))
        mastername = listentry%field%name

        call list_index (fwrtpr(1,t), mastername, ff)
        if (ff > 0) then
          prec_wrt = getflag(fwrtpr(ff,t))
        else
          prec_wrt = ' '
        end if

        call list_index (fincl(1,t), mastername, ff)

        if (ff > 0) then
          avgflag = getflag (fincl(ff,t))
          call inifld (t, listentry, avgflag,  prec_wrt)
        else if ((.not. empty_htapes) .or. (is_initfile(file_index=t))) then
          call list_index (fexcl(1,t), mastername, ff)
          if (ff == 0 .and. listentry%actflag(t)) then
            call inifld (t, listentry, ' ', prec_wrt)
          else
            listentry%actflag(t) = .false.
          end if
        else
          listentry%actflag(t) = .false.
        end if
        listentry=>listentry%next_entry

      end do
      !
      ! If column output is specified make sure there are some fields defined
      ! for that tape
      !
      if (nflds(t) .eq. 0 .and. fincllonlat(1,t) .ne. ' ') then
        write(errormsg,'(a,i2,a)') 'FLDLST: Column output is specified for tape ',t,' but no fields defined for that tape.'
        call endrun(errormsg)
      else
        call patch_init(t)
      end if
      !
      ! Specification of tape contents now complete.  Sort each list of active
      ! entries for efficiency in OUTFLD.  Simple bubble sort.
      !
!!XXgoldyXX: v In the future, we will sort according to decomp to speed I/O
      do f=nflds(t)-1,1,-1
        do ff=1,f

          if (tape(t)%hlist(ff)%field%name > tape(t)%hlist(ff+1)%field%name) then

            tmp = tape(t)%hlist(ff)
            tape(t)%hlist(ff  ) = tape(t)%hlist(ff+1)
            tape(t)%hlist(ff+1) = tmp

          else if (tape(t)%hlist(ff  )%field%name == tape(t)%hlist(ff+1)%field%name) then

            write(errormsg,'(2a,2(a,i3))') 'FLDLST: Duplicate field: ', &
                 trim(tape(t)%hlist(ff)%field%name),', tape = ', t, ', ff = ', ff
            call endrun(errormsg)

          end if

        end do
      end do

    end do    ! do t=1,ptapes
    deallocate(gridsontape)

    call print_active_fldlst()

    !
    ! Packing density, ndens: With netcdf, only 1 (nf_double) and 2 (pio_real)
    ! are allowed
    !
    do t=1,ptapes
      if (ndens(t) == 1) then
        ncprec(t) = pio_double
      else if (ndens(t) == 2) then
        ncprec(t) = pio_real
      else
        call endrun ('FLDLST: ndens must be 1 or 2')
      end if

    end do
    !
    !  Now that masterlinkedlist is defined, construct primary and secondary hashing
    !  tables.
    !
    call bld_outfld_hash_tbls()
    call bld_htapefld_indices()

    return
  end subroutine fldlst

!#########################################################################################

subroutine print_active_fldlst()

   integer :: f, ff, i, t
   integer :: num_patches

   character(len=6) :: prec_str
   character(len=max_chars) :: fldname, fname_tmp

   type(active_entry), pointer :: hfile(:) => null()  ! history files

   if (masterproc) then

      hfile=>history_tape

      do t=1,ptapes

         if (nflds(t) > 0) then
            write(iulog,*) ' '
            write(iulog,*)'FLDLST: History file ', t, ' contains ', nflds(t), ' fields'

            if (is_initfile(file_index=t)) then
               write(iulog,*) ' Write frequency:                 ',inithist,' (INITIAL CONDITIONS)'
            else
               if (nhtfrq(t) == 0) then
                  write(iulog,*) ' Write frequency:                  MONTHLY'
               else
                  write(iulog,*) ' Write frequency:                 ',nhtfrq(t)
               end if
            end if

            write(iulog,*) ' Filename specifier:              ', trim(hfilename_spec(t))

            prec_str = 'double'
            if (ndens(t) == 2) prec_str = 'single'
            write(iulog,*) ' Output precision:                ', prec_str
            write(iulog,*) ' Number of time samples per file: ', mfilt(t)

            ! grid info
            if (associated(hfile(t)%patches)) then
               write(iulog,*) ' Fields are represented on columns (FIELD_LON_LAT)'
            else if (associated(hfile(t)%grid_ids)) then
               write(iulog,*) ' Fields are represented on global grids:'
               do i = 1, size(hfile(t)%grid_ids)
                  write(iulog,*) ' ', hfile(t)%grid_ids(i)
               end do
            else
               call endrun('print_active_fldlst: error in active_entry object')
            end if

            write(iulog,*)' Included fields are:'

         end if

         do f = 1, nflds(t)
            if (associated(hfile(t)%patches)) then
               num_patches = size(hfile(t)%patches)
               fldname = strip_suffix(hfile(t)%hlist(f)%field%name)
               do i = 1, num_patches
                  ff = (f-1)*num_patches + i
                  fname_tmp = trim(fldname)
                  call hfile(t)%patches(i)%field_name(fname_tmp)
                  write(iulog,9000) ff, fname_tmp, hfile(t)%hlist(f)%field%units, &
                     hfile(t)%hlist(f)%field%numlev, hfile(t)%hlist(f)%avgflag,   &
                     trim(hfile(t)%hlist(f)%field%long_name)
               end do
            else
               fldname = hfile(t)%hlist(f)%field%name
               write(iulog,9000) f, fldname, hfile(t)%hlist(f)%field%units,  &
                  hfile(t)%hlist(f)%field%numlev, hfile(t)%hlist(f)%avgflag, &
                  trim(hfile(t)%hlist(f)%field%long_name)
            end if

         end do

      end do

   end if

9000 format(i5, 1x, a32, 1x, a16, 1x, i4, 1x, a1, 2x, 256a)

end subroutine print_active_fldlst

!#########################################################################################

  subroutine inifld (t, listentry, avgflag, prec_wrt)
    use cam_grid_support, only: cam_grid_is_zonal
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Add a field to the active list for a history tape
    !
    ! Method: Copy the data from the master field list to the active list for the tape
    !         Also: define mapping arrays from (col,chunk) -> (lon,lat)
    !
    ! Author: CCM Core Group
    !
    !-----------------------------------------------------------------------


    !
    ! Arguments
    !
    integer, intent(in)          :: t         ! history tape index

    type(master_entry), pointer  :: listentry

    character(len=1), intent(in) :: avgflag   ! averaging flag
    character(len=1), intent(in) :: prec_wrt  ! history output precision flag
    !
    ! Local workspace
    !
    integer :: n                  ! field index on defined tape


    !
    ! Ensure that it is not to late to add a field to the history tape
    !
    if (htapes_defined) then
      call endrun ('INIFLD: Attempt to add field '//listentry%field%name//' after history files set')
    end if

    nflds(t) = nflds(t) + 1
    n = nflds(t)
    !
    ! Copy field info.
    !
    if(n > size(tape(t)%hlist)) then
      write(iulog,*) 'tape field miscount error ', n, size(tape(t)%hlist)
      call endrun()
    end if

    tape(t)%hlist(n)%field = listentry%field

    select case (prec_wrt)
    case (' ')
      if (ndens(t) == 1) then
        tape(t)%hlist(n)%hwrt_prec = 8
      else
        tape(t)%hlist(n)%hwrt_prec = 4
      end if
    case ('4')
      tape(t)%hlist(n)%hwrt_prec = 4
      if (masterproc) then
        write(iulog,*) 'INIFLD: Output data type for ', tape(t)%hlist(n)%field%name, &
             ' is real*4'
      end if
    case ('8')
      tape(t)%hlist(n)%hwrt_prec = 8
      if (masterproc) then
        write(iulog,*) 'INIFLD: Output data type for ', tape(t)%hlist(n)%field%name, &
             ' is real*8'
      end if
    case default
      call endrun ('INIFLD: unknown prec_wrt='//prec_wrt)
    end select
    !
    ! Override the default averaging (masterlist) averaging flag if non-blank
    !
    if (avgflag == ' ') then
      tape(t)%hlist(n)%avgflag = listentry%avgflag(t)
      tape(t)%hlist(n)%time_op = listentry%time_op(t)
    else
      tape(t)%hlist(n)%avgflag = avgflag
      call AvgflagToString(avgflag, tape(t)%hlist(n)%time_op)
    end if

    ! Some things can't be done with zonal fields
    if (cam_grid_is_zonal(listentry%field%decomp_type)) then
      if (tape(t)%hlist(n)%avgflag == 'L') then
        call endrun("Cannot perform local time processing on zonal data ("//trim(listentry%field%name)//")")
      else if (is_satfile(t)) then
        call endrun("Zonal data not valid for satellite history ("//trim(listentry%field%name)//")")
      end if
    end if

#ifdef HDEBUG
    if (masterproc) then
      write(iulog,'(a,i0,3a,i0,a,i2)')'HDEBUG: ',__LINE__,' field ',          &
           trim(tape(t)%hlist(n)%field%name), ' added as field number ', n,   &
           ' on tape ', t
      write(iulog,'(2a)')'  units     = ',trim(tape(t)%hlist(n)%field%units)
      write(iulog,'(a,i0)')'  numlev    = ',tape(t)%hlist(n)%field%numlev
      write(iulog,'(2a)')'  avgflag   = ',tape(t)%hlist(n)%avgflag
      write(iulog,'(3a)')'  time_op   = "',trim(tape(t)%hlist(n)%time_op),'"'
      write(iulog,'(a,i0)')'  hwrt_prec = ',tape(t)%hlist(n)%hwrt_prec
    end if
#endif

    return
  end subroutine inifld


  subroutine patch_init(t)
    use cam_history_support, only: history_patch_t
    use cam_grid_support,    only: cam_grid_compute_patch

    ! Dummy arguments
    integer, intent(in)               :: t     ! Current tape

    ! Local variables
    integer                           :: ff    ! Loop over fincllonlat entries
    integer                           :: i     ! General loop index
    integer                           :: npatches
    type(history_patch_t), pointer    :: patchptr

    character(len=max_chars)          :: errormsg
    character(len=max_chars)          :: lonlatname(pflds)
    real(r8)                          :: beglon, beglat, endlon, endlat

    !
    ! Setup column information if this field will be written as group
    ! First verify the column information in the namelist
    ! Duplicates are an error, but we can just ignore them
    !

    ! I know, this shouldn't happen . . . yet: (better safe than sorry)
    if (associated(tape(t)%patches)) then
      do i = 1, size(tape(t)%patches)
        call tape(t)%patches(i)%deallocate()
      end do
      deallocate(tape(t)%patches)
      nullify(tape(t)%patches)
    end if

    ! First, count the number of patches and check for duplicates
    ff = 1  ! Index of fincllonlat entry
    npatches = 0   ! Number of unique patches in namelist entry
    do while (len_trim(fincllonlat(ff, t)) > 0)
      npatches = npatches + 1
      lonlatname(npatches) = trim(fincllonlat(ff, t))
      ! Check for duplicates
      do i = 1, npatches - 1
        if (trim(lonlatname(i)) == trim(lonlatname(npatches))) then
          write(errormsg, '(a,i0,3a)') 'Duplicate fincl', t, 'lonlat entry.', &
               'Duplicate entry is ', trim(lonlatname(i))
          write(iulog, *) 'patch_init: WARNING: '//errormsg
          ! Remove the new entry
          lonlatname(npatches) = ''
          npatches = npatches - 1
          exit
        end if
      end do
      ff = ff + 1
    end do

    ! Now we know how many patches, allocate space
    if (npatches > 0) then
      if (collect_column_output(t)) then
        allocate(tape(t)%patches(1))
      else
        allocate(tape(t)%patches(npatches))
      end if

      ! For each lat/lon specification, parse and create a patch for each grid
      do ff = 1, npatches
        if (collect_column_output(t)) then
          ! For colleccted column output, we only have one patch
          patchptr => tape(t)%patches(1)
        else
          patchptr => tape(t)%patches(ff)
          patchptr%namelist_entry = trim(lonlatname(ff))
        end if
        ! We need to set up one patch per (active) grid
        patchptr%collected_output = collect_column_output(t)
        call parseLonLat(lonlatname(ff),                                      &
             beglon, endlon, patchptr%lon_axis_name,                          &
             beglat, endlat, patchptr%lat_axis_name)
        if (associated(patchptr%patches)) then
          ! One last sanity check
          if (.not. collect_column_output(t)) then
            write(errormsg, '(a,i0,2a)') 'Attempt to overwrite fincl', t,     &
                 'lonlat entry, ', trim(patchptr%namelist_entry)
            call endrun('patch_init: '//errormsg)
          end if
        else
          allocate(patchptr%patches(size(tape(t)%grid_ids)))
        end if
        do i = 1, size(tape(t)%grid_ids)
          call cam_grid_compute_patch(tape(t)%grid_ids(i), patchptr%patches(i),&
               beglon, endlon, beglat, endlat, collect_column_output(t))
        end do
        nullify(patchptr)
      end do
    end if
    ! We are done processing this tape's fincl#lonlat entries. Now,
    ! compact each patch so that the output variables have no holes
    ! We wait until now for when collect_column_output(t) is .true. since
    !    all the fincl#lonlat entries are concatenated
    if (associated(tape(t)%patches)) then
      do ff = 1, size(tape(t)%patches)
        call tape(t)%patches(ff)%compact()
      end do
    end if

  end subroutine patch_init

  !#######################################################################

  subroutine strip_null(str)
    character(len=*), intent(inout) :: str
    integer :: i
    do i=1,len(str)
      if(ichar(str(i:i))==0) str(i:i)=' '
    end do
  end subroutine strip_null

  character(len=max_fieldname_len) function strip_suffix (name)
    !
    !----------------------------------------------------------
    !
    ! Purpose:  Strip "&IC" suffix from fieldnames if it exists
    !
    !----------------------------------------------------------
    !
    ! Arguments
    !
    character(len=*), intent(in) :: name
    !
    ! Local workspace
    !
    integer :: n
    !
    !-----------------------------------------------------------------------
    !
    strip_suffix = ' '

    do n = 1,fieldname_len
      strip_suffix(n:n) = name(n:n)
      if(name(n+1:n+1         ) == ' '                       ) return
      if(name(n+1:n+fieldname_suffix_len) == fieldname_suffix) return
    end do

    strip_suffix(fieldname_len+1:max_fieldname_len) = name(fieldname_len+1:max_fieldname_len)

    return

  end function strip_suffix

  !#######################################################################

  character(len=fieldname_len) function getname (inname)
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: retrieve name portion of inname
    !
    ! Method:  If an averaging flag separater character is present (":") in inname,
    !          lop it off
    !
    !-------------------------------------------------------------------------------
    !
    ! Arguments
    !
    character(len=*), intent(in) :: inname
    !
    ! Local workspace
    !
    integer :: length
    integer :: i

    length = len (inname)

    if (length < fieldname_len .or. length > fieldname_lenp2) then
      write(iulog,*) 'GETNAME: bad length=',length
      call endrun
    end if

    getname = ' '
    do i=1,fieldname_len
      if (inname(i:i) == ':') exit
      getname(i:i) = inname(i:i)
    end do

    return
  end function getname

  !#######################################################################

  ! parseRangeString: Parse either a coordinate descriptor (e.g., 10S) or a
  !                   coordinate range (e.g., 10e:20e)
  !                   chars represents the allowed coordinate character.
  !                   NB: Does not validate numerical values (e.g., lat <= 90)
  subroutine parseRangeString(rangestr, chars, begval, begchar, begname, endval, endchar, endname)

    ! Dummy arguments
    character(len=*),       intent(in)    :: rangestr
    character(len=*),       intent(in)    :: chars
    real(r8),               intent(out)   :: begval
    character,              intent(out)   :: begchar
    character(len=*),       intent(out)   :: begname
    real(r8),               intent(out)   :: endval
    character,              intent(out)   :: endchar
    character(len=*),       intent(out)   :: endname

    ! Local variables
    character(len=128)                    :: errormsg
    integer                               :: colonpos
    integer                               :: beglen, endlen

    ! First, see if we have a position or a range
    colonpos = scan(rangestr, ':')
    if (colonpos == 0) then
      begname = trim(rangestr)
      beglen = len_trim(begname)
      endname = trim(begname)
    else
      beglen = colonpos - 1
      begname = rangestr(1:beglen)
      endname = trim(rangestr(colonpos+1:))
      endlen = len_trim(endname)
    end if
    ! begname should be a number (integer or real) followed by a character
    if (verify(begname, '0123456789.') /= beglen) then
      write(errormsg, *) 'Coordinate range must begin with number, ', begname
      call endrun('parseRangeString: '//errormsg)
    end if
    if (verify(begname(beglen:beglen), chars) /= 0) then
      write(errormsg, *) 'Coordinate range must end with character in the ',  &
           'set [', trim(chars), '] ', begname
      call endrun('parseRangeString: '//errormsg)
    end if
    ! begname parses so collect the values
    read(begname(1:beglen-1), *) begval
    begchar = begname(beglen:beglen)
    if (colonpos /= 0) then
      ! endname should be a number (integer or real) followed by a character
      if (verify(endname, '0123456789.') /= endlen) then
        write(errormsg, *) 'Coordinate range must begin with number, ', endname
        call endrun('parseRangeString: '//errormsg)
      end if
      if (verify(endname(endlen:endlen), chars) /= 0) then
        write(errormsg, *) 'Coordinate range must end with character in the ',&
             'set [', trim(chars), '] ', endname
        call endrun('parseRangeString: '//errormsg)
      end if
      ! endname parses so collect the values
      read(endname(1:endlen-1), *) endval
      endchar = endname(endlen:endlen)
    else
      endval = begval
      endchar = begchar
    end if

  end subroutine parseRangeString

  ! parseLonLat: Parse a lon_lat description allowed by the fincllonlat(n)
  !              namelist entries. Returns the starting and ending values of
  !              the point or range specified.
  !              NB: Does not validate the range against any particular grid
  subroutine parseLonLat(lonlatname, beglon, endlon, lonname, beglat, endlat, latname)

    ! Dummy arguments
    character(len=*),       intent(in)    :: lonlatname
    real(r8),               intent(out)   :: beglon
    real(r8),               intent(out)   :: endlon
    character(len=*),       intent(out)   :: lonname
    real(r8),               intent(out)   :: beglat
    real(r8),               intent(out)   :: endlat
    character(len=*),       intent(out)   :: latname

    ! Local variables
    character(len=128)                    :: errormsg
    character(len=MAX_CHARS)              :: lonstr, latstr
    character(len=MAX_CHARS)              :: begname, endname
    character                             :: begchar, endchar
    integer                               :: underpos

    !
    ! make sure _ separator is present
    !
    underpos = scan(lonlatname, '_')
    if (underpos == 0) then
      write(errormsg,*) 'Improperly formatted fincllonlat string. ',          &
           'Missing underscore character (xxxE_yyyS) ', lonlatname
      call endrun('parseLonLat: '//errormsg)
    end if

    ! Break out the longitude and latitude sections
    lonstr = lonlatname(:underpos-1)
    latstr = trim(lonlatname(underpos+1:))

    ! Parse the longitude section
    call parseRangeString(lonstr, 'eEwW', beglon, begchar, begname, endlon, endchar, endname)
    ! Convert longitude to degrees East
    if ((begchar == 'w') .or. (begchar == 'W')) then
      if (beglon > 0.0_r8) then
        beglon = 360._r8 - beglon
      end if
    end if
    if ((beglon < 0._r8) .or. (beglon > 360._r8)) then
      write(errormsg, *) 'Longitude specification out of range, ', trim(begname)
      call endrun('parseLonLat: '//errormsg)
    end if
    if ((endchar == 'w') .or. (endchar == 'W')) then
      if (endlon > 0.0_r8) then
        endlon = 360._r8 - endlon
      end if
    end if
    if ((endlon < 0._r8) .or. (endlon > 360._r8)) then
      write(errormsg, *) 'Longitude specification out of range, ', trim(endname)
      call endrun('parseLonLat: '//errormsg)
    end if
    if (beglon == endlon) then
      lonname = trim(begname)
    else
      lonname = trim(begname)//'_to_'//trim(endname)
    end if

    ! Parse the latitude section
    call parseRangeString(latstr, 'nNsS', beglat, begchar, begname, endlat, endchar, endname)
    ! Convert longitude to degrees East
    if ((begchar == 's') .or. (begchar == 'S')) then
      beglat = (-1._r8) * beglat
    end if
    if ((beglat < -90._r8) .or. (beglat > 90._r8)) then
      write(errormsg, *) 'Latitude specification out of range, ', trim(begname)
      call endrun('parseLonLat: '//errormsg)
    end if
    if ((endchar == 's') .or. (endchar == 'S')) then
      endlat = (-1._r8) * endlat
    end if
    if ((endlat < -90._r8) .or. (endlat > 90._r8)) then
      write(errormsg, *) 'Latitude specification out of range, ', trim(endname)
      call endrun('parseLonLat: '//errormsg)
    end if
    if (beglat == endlat) then
      latname = trim(begname)
    else
      latname = trim(begname)//'_to_'//trim(endname)
    end if

  end subroutine parseLonLat


  !#######################################################################

  character(len=1) function getflag (inname)
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: retrieve flag portion of inname
    !
    ! Method:  If an averaging flag separater character is present (":") in inname,
    !          return the character after it as the flag
    !
    !-------------------------------------------------------------------------------
    !
    ! Arguments
    !
    character(len=*), intent(in) :: inname   ! character string
    !
    ! Local workspace
    !
    integer :: length         ! length of inname
    integer :: i              ! loop index

    length = len (inname)

    if (length /= fieldname_lenp2) then
      write(iulog,*) 'GETFLAG: bad length=',length
      call endrun
    end if

    getflag = ' '
    do i=1,fieldname_lenp2-1
      if (inname(i:i) == ':') then
        getflag = inname(i+1:i+1)
        exit
      end if
    end do

    return
  end function getflag

  !#######################################################################

  subroutine list_index (list, name, index)
    !
    ! Input arguments
    !
    character(len=*), intent(in) :: list(pflds) ! input list of names, possibly ":" delimited
    character(len=max_fieldname_len), intent(in) :: name ! name to be searched for
    !
    ! Output arguments
    !
    integer, intent(out) :: index               ! index of "name" in "list"
    !
    ! Local workspace
    !
    character(len=fieldname_len) :: listname    ! input name with ":" stripped off.
    integer f                       ! field index

    index = 0
    do f=1,pflds
      !
      ! Only list items
      !
      listname = getname (list(f))
      if (listname == ' ') exit
      if (listname == name) then
        index = f
        exit
      end if
    end do

    return
  end subroutine list_index

  !#######################################################################

  recursive subroutine outfld (fname, field, idim, c, avg_subcol_field)
    use cam_history_buffers, only: hbuf_accum_inst, hbuf_accum_add, hbuf_accum_variance,  &
         hbuf_accum_add00z, hbuf_accum_max, hbuf_accum_min,          &
         hbuf_accum_addlcltime
    use cam_history_support, only: dim_index_2d
    use subcol_pack_mod,     only: subcol_unpack
    use cam_grid_support,    only: cam_grid_id

    interface
      subroutine subcol_field_avg_handler(idim, field_in, c, field_out)
        use shr_kind_mod, only: r8 => shr_kind_r8
        integer,  intent(in)  :: idim
        real(r8), intent(in)  :: field_in(idim, *)
        integer,  intent(in)  :: c
        real(r8), intent(out) :: field_out(:,:)
      end subroutine subcol_field_avg_handler
    end interface

    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Accumulate (or take min, max, etc. as appropriate) input field
    !          into its history buffer for appropriate tapes
    !
    ! Method: Check 'masterlist' whether the requested field 'fname' is active
    !         on one or more history tapes, and if so do the accumulation.
    !         If not found, return silently.
    !         subcol_field_avg_handler:
    !            An interface into subcol_field_avg without creating a dependency as
    !            this would cause a dependency loop. See subcol.F90
    ! Note: We cannot know a priori if field is a grid average field or a subcolumn
    !       field because many fields passed to outfld are defined on ncol rather
    !       than pcols or psetcols. Therefore, we use the avg_subcol_field input
    !       to determine whether to average the field input before accumulation.
    !       NB: If output is on a subcolumn grid (requested in addfle), it is
    !           an error to use avg_subcol_field. A subcolumn field is assumed and
    !           subcol_unpack is called before accumulation.
    !
    ! Author: CCM Core Group
    !
    !-----------------------------------------------------------------------
    !
    ! Arguments
    !
    character(len=*), intent(in) :: fname ! Field name--should be 8 chars long

    ! For structured grids, idim is the local longitude dimension.
    ! For unstructured grids, idim is the local column dimension
    ! For phys_decomp, it should be pcols or pcols*psubcols
    integer, intent(in)           :: idim
    real(r8), intent(in)          :: field(idim,*) ! Array containing field values
    integer, intent(in)           :: c             ! chunk (physics) or latitude (dynamics) index
    logical, optional, intent(in) :: avg_subcol_field
    !
    ! Local variables
    !
    integer               :: t, f          ! tape, field indices

    character*1           :: avgflag       ! averaging flag

    type (active_entry), pointer :: otape(:) ! Local history_tape pointer
    real(r8),pointer      :: hbuf(:,:)     ! history buffer
    real(r8),pointer      :: sbuf(:,:)     ! variance buffer
    integer, pointer      :: nacs(:)       ! accumulation counter
    integer               :: begdim2, enddim2, endi
    integer               :: phys_decomp
    type (dim_index_2d)   :: dimind        ! 2-D dimension index
    logical               :: flag_xyfill   ! non-applicable xy points flagged with fillvalue
    real(r8)              :: fillvalue
    real(r8), allocatable :: afield(:,:)   ! Averaged field values
    real(r8), allocatable :: ufield(:,:,:) ! Unpacked field values
    integer               :: ff            ! masterlist index pointer
    integer               :: i, j
    logical               :: found
    logical               :: avg_subcols   ! average subcols before accum
    !-----------------------------------------------------------------------

    call get_field_properties(fname, found, tape_out=otape, ff_out=ff)
    phys_decomp = cam_grid_id('physgrid')

    ! If this field is not active, return now
    if (.not. found) then
      return
    end if

    !
    ! Note, the field may be on any or all of the history files (primary
    ! and auxiliary).
    !
    !      write(iulog,*)'fname_loc=',fname_loc
    do t = 1, ptapes
      if ( .not. masterlist(ff)%thisentry%actflag(t)) cycle
      f = masterlist(ff)%thisentry%htapeindx(t)
      !
      ! Update history buffer
      !
      flag_xyfill = otape(t)%hlist(f)%field%flag_xyfill
      fillvalue = otape(t)%hlist(f)%field%fillvalue
      avgflag = otape(t)%hlist(f)%avgflag
      nacs   => otape(t)%hlist(f)%nacs(:,c)
      hbuf => otape(t)%hlist(f)%hbuf(:,:,c)
      if (associated(tape(t)%hlist(f)%sbuf)) then
         sbuf => otape(t)%hlist(f)%sbuf(:,:,c)
      endif
      dimind = otape(t)%hlist(f)%field%get_dims(c)

      ! See notes above about validity of avg_subcol_field
      if (otape(t)%hlist(f)%field%is_subcol) then
        if (present(avg_subcol_field)) then
          call endrun('OUTFLD: Cannot average '//trim(fname)//', subcolumn output was requested in addfld')
        end if
        avg_subcols = .false.
      else if (otape(t)%hlist(f)%field%decomp_type == phys_decomp) then
        if (present(avg_subcol_field)) then
          avg_subcols = avg_subcol_field
        else
          avg_subcols = .false.
        end if
      else ! Any dynamics decomposition
        if (present(avg_subcol_field)) then
          call endrun('OUTFLD: avg_subcol_field only valid for physgrid')
        else
          avg_subcols = .false.
        end if
      end if

      begdim2 = otape(t)%hlist(f)%field%begdim2
      enddim2 = otape(t)%hlist(f)%field%enddim2
      if (avg_subcols) then
        allocate(afield(pcols, begdim2:enddim2))
        call subcol_field_avg_handler(idim, field, c, afield)
        ! Hack! Avoid duplicating select statement below
        call outfld(fname, afield, pcols, c)
        deallocate(afield)
      else if (otape(t)%hlist(f)%field%is_subcol) then
        ! We have to assume that using mdimnames (e.g., psubcols) is
        ! incompatible with the begdimx, enddimx usage (checked in addfld)
        ! Since psubcols is included in levels, take that out
        endi = (enddim2 - begdim2 + 1) / psubcols
        allocate(ufield(pcols, psubcols, endi))
        allocate(afield(pcols*psubcols, endi))
        do j = 1, endi
          do i = 1, idim
            afield(i, j) = field(i, j)
          end do
        end do
        ! Initialize unused aray locations.
        if (idim < pcols*psubcols) then
          if (flag_xyfill) then
            afield(idim+1:pcols*psubcols, :) = fillvalue
          else
            afield(idim+1:pcols*psubcols, :) = 0.0_r8
          end if
        end if
        if (flag_xyfill) then
          call subcol_unpack(c, afield, ufield, fillvalue)
        else
          call subcol_unpack(c, afield, ufield)
        end if
        deallocate(afield)
        select case (avgflag)

        case ('I') ! Instantaneous
          call hbuf_accum_inst(hbuf, ufield, nacs, dimind, pcols,        &
               flag_xyfill, fillvalue)

        case ('A') ! Time average
          call hbuf_accum_add(hbuf, ufield, nacs, dimind, pcols,         &
               flag_xyfill, fillvalue)

        case ('B') ! Time average only 00z values
          call hbuf_accum_add00z(hbuf, ufield, nacs, dimind, pcols,      &
               flag_xyfill, fillvalue)

        case ('X') ! Maximum over time
          call hbuf_accum_max (hbuf, ufield, nacs, dimind, pcols,        &
               flag_xyfill, fillvalue)

        case ('M') ! Minimum over time
          call hbuf_accum_min(hbuf, ufield, nacs, dimind, pcols,         &
               flag_xyfill, fillvalue)

        case ('L')
          call hbuf_accum_addlcltime(hbuf, ufield, nacs, dimind, pcols,  &
               flag_xyfill, fillvalue, c,                                &
               otape(t)%hlist(f)%field%decomp_type,                      &
               lcltod_start(t), lcltod_stop(t))

        case ('S') ! Standard deviation
          call hbuf_accum_variance(hbuf, sbuf, ufield, nacs, dimind, pcols,&
               flag_xyfill, fillvalue)

        case default
          call endrun ('OUTFLD: invalid avgflag='//avgflag)

        end select
        deallocate(ufield)
      else
        select case (avgflag)

        case ('I') ! Instantaneous
          call hbuf_accum_inst(hbuf, field, nacs, dimind, idim,          &
               flag_xyfill, fillvalue)

        case ('A') ! Time average
          call hbuf_accum_add(hbuf, field, nacs, dimind, idim,           &
               flag_xyfill, fillvalue)

        case ('B') ! Time average only 00z values
          call hbuf_accum_add00z(hbuf, field, nacs, dimind, idim,        &
               flag_xyfill, fillvalue)

        case ('X') ! Maximum over time
          call hbuf_accum_max (hbuf, field, nacs, dimind, idim,          &
               flag_xyfill, fillvalue)

        case ('M') ! Minimum over time
          call hbuf_accum_min(hbuf, field, nacs, dimind, idim,           &
               flag_xyfill, fillvalue)

        case ('L')
          call hbuf_accum_addlcltime(hbuf, field, nacs, dimind, idim,    &
               flag_xyfill, fillvalue, c,                                &
               otape(t)%hlist(f)%field%decomp_type,                      &
               lcltod_start(t), lcltod_stop(t))

        case ('S') ! Standard deviation
          call hbuf_accum_variance(hbuf, sbuf, field, nacs, dimind, idim,&
               flag_xyfill, fillvalue)

        case default
          call endrun ('OUTFLD: invalid avgflag='//avgflag)

        end select
      end if

    end do

    return
  end subroutine outfld

  !#######################################################################

  subroutine get_field_properties(fname, found, tape_out, ff_out, no_tape_check_in)

    implicit none
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: If fname is active, lookup and return field information
    !
    ! Method: Check 'masterlist' whether the requested field 'fname' is active
    !         on one or more history tapes, and if so, return the requested
    !         field information
    !
    ! Author: goldy
    !
    !-----------------------------------------------------------------------
    !
    ! Arguments
    !
    character(len=*),   intent(in)  :: fname ! Field name--should be 8 chars long
    logical,            intent(out) :: found ! Set to true if fname is active
    type(active_entry), pointer, optional :: tape_out(:)
    integer,            intent(out), optional :: ff_out
    logical,            intent(in), optional  :: no_tape_check_in

    !
    ! Local variables
    !
    character*(max_fieldname_len) :: fname_loc  ! max-char equivalent of fname
    integer :: t, ff          ! tape, masterindex indices
    logical :: no_tape_check
     !-----------------------------------------------------------------------

    ! Need to re-cast the field name so that the hashing works #hackalert
    fname_loc = fname
    ff = get_masterlist_indx(fname_loc)

    ! Set the no_tape_check to false, unless is passed in
    if (present(no_tape_check_in)) then
       no_tape_check = no_tape_check_in
    else
       no_tape_check = .false.
    end if

    ! Set found to .false. so we can return early if fname is not active
    found = .false.
    if (present(tape_out)) then
      nullify(tape_out)
    end if
    if (present(ff_out)) then
      ff_out = -1
    end if

    !
    !  If ( ff < 0 ), the field is not defined on the masterlist. This check
    !  is necessary because of coding errors calling outfld without first defining
    !  the field on masterlist.
    !
    if ( ff < 0 ) then
      return
    end if
    !
    !  Next, check to see whether this field is active on one or more history
    !  tapes.
    !
    if (no_tape_check) then
      if (present(ff_out)) ff_out   =  ff  ! Set the output index and return without checking tapes
      return
    else if ( .not. masterlist(ff)%thisentry%act_sometape )  then
      return
    end if
    !
    ! Note, the field may be on any or all of the history files (primary
    ! and auxiliary).
    !

    do t=1, ptapes
      if (masterlist(ff)%thisentry%actflag(t)) then
        found    =  .true.
        if (present(tape_out)) then
          tape_out => history_tape
        end if
        if (present(ff_out)) then
          ff_out   =  ff
        end if
        ! We found the info so we are done with the loop
        exit
      end if
    end do

  end subroutine get_field_properties

  !#######################################################################

  logical function is_initfile (file_index)
    !
    !------------------------------------------------------------------------
    !
    ! Purpose: to determine:
    !
    !   a) if an IC file is active in this model run at all
    !       OR,
    !   b) if it is active, is the current file index referencing the IC file
    !      (IC file is always at ptapes)
    !
    !------------------------------------------------------------------------
    !
    ! Arguments
    !
    integer, intent(in), optional :: file_index ! index of file in question

    is_initfile = .false.

    if (present(file_index)) then
      if (inithist /= 'NONE' .and. file_index == ptapes) is_initfile = .true.
    else
      if (inithist /= 'NONE'                           ) is_initfile = .true.
    end if

    return

  end function is_initfile

  !#######################################################################

  integer function strcmpf (name1, name2)
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Return the lexical difference between two strings
    !
    ! Method: Use ichar() intrinsic as we loop through the names
    !
    !-----------------------------------------------------------------------
    !
    ! Arguments
    !
    character(len=max_fieldname_len), intent(in) :: name1, name2 ! strings to compare
    integer n                                     ! loop index

    do n=1,max_fieldname_len
      strcmpf = ichar(name1(n:n)) - ichar(name2(n:n))
      if (strcmpf /= 0) exit
    end do

    return
  end function strcmpf

  !#######################################################################

  subroutine h_inquire (t)
    use pio,           only: pio_inq_varid, pio_inq_attlen
    use cam_pio_utils, only: cam_pio_handle_error
   !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Ensure that the proper variables are on a history file
    !
    ! Method: Issue the appropriate netcdf wrapper calls
    !
    !-----------------------------------------------------------------------
    !
    ! Arguments
    !
    integer, intent(in) :: t   ! tape index
    !
    ! Local workspace
    !
    integer                  :: f            ! field index
    integer                  :: ierr
    integer                  :: i
    integer                  :: num_patches
    integer(pio_offset_kind) :: mdimsize
    character(len=max_chars) :: fldname, fname_tmp, basename

    !
    !
    ! Dimension id's
    !
    tape => history_tape



    !
    ! Create variables for model timing and header information
    !
    if(.not. is_satfile(t)) then
      ierr=pio_inq_varid (tape(t)%File,'ndcur   ',    tape(t)%ndcurid)
      ierr=pio_inq_varid (tape(t)%File,'nscur   ',    tape(t)%nscurid)
      ierr=pio_inq_varid (tape(t)%File,'nsteph  ',    tape(t)%nstephid)

      ierr=pio_inq_varid (tape(t)%File,'time_bnds',   tape(t)%tbndid)
      ierr=pio_inq_varid (tape(t)%File,'date_written',tape(t)%date_writtenid)
      ierr=pio_inq_varid (tape(t)%File,'time_written',tape(t)%time_writtenid)
#if ( defined BFB_CAM_SCAM_IOP )
      ierr=pio_inq_varid (tape(t)%File,'tsec    ',tape(t)%tsecid)
      ierr=pio_inq_varid (tape(t)%File,'bdate   ',tape(t)%bdateid)
#endif
      if (.not. is_initfile(file_index=t) ) then
        ! Don't write the GHG/Solar forcing data to the IC file.  It is never
        ! read from that file so it's confusing to have it there.
        ierr=pio_inq_varid (tape(t)%File,'co2vmr  ',    tape(t)%co2vmrid)
        ierr=pio_inq_varid (tape(t)%File,'ch4vmr  ',    tape(t)%ch4vmrid)
        ierr=pio_inq_varid (tape(t)%File,'n2ovmr  ',    tape(t)%n2ovmrid)
        ierr=pio_inq_varid (tape(t)%File,'f11vmr  ',    tape(t)%f11vmrid)
        ierr=pio_inq_varid (tape(t)%File,'f12vmr  ',    tape(t)%f12vmrid)
        ierr=pio_inq_varid (tape(t)%File,'sol_tsi ',    tape(t)%sol_tsiid)
        if (solar_parms_on) then
          ierr=pio_inq_varid (tape(t)%File,'f107    ',    tape(t)%f107id)
          ierr=pio_inq_varid (tape(t)%File,'f107a   ',    tape(t)%f107aid)
          ierr=pio_inq_varid (tape(t)%File,'f107p   ',    tape(t)%f107pid)
          ierr=pio_inq_varid (tape(t)%File,'kp      ',    tape(t)%kpid)
          ierr=pio_inq_varid (tape(t)%File,'ap      ',    tape(t)%apid)
        endif
        if (solar_wind_on) then
          ierr=pio_inq_varid (tape(t)%File,'byimf', tape(t)%byimfid)
          ierr=pio_inq_varid (tape(t)%File,'bzimf', tape(t)%bzimfid)
          ierr=pio_inq_varid (tape(t)%File,'swvel', tape(t)%swvelid)
          ierr=pio_inq_varid (tape(t)%File,'swden', tape(t)%swdenid)
        endif
      end if
    end if
    ierr=pio_inq_varid (tape(t)%File,'date    ',    tape(t)%dateid)
    ierr=pio_inq_varid (tape(t)%File,'datesec ',    tape(t)%datesecid)
    ierr=pio_inq_varid (tape(t)%File,'time    ',    tape(t)%timeid)


    !
    ! Obtain variable name from ID which was read from restart file
    !
    do f=1,nflds(t)
      if(.not. associated(tape(t)%hlist(f)%varid)) then
        if (associated(tape(t)%patches)) then
          allocate(tape(t)%hlist(f)%varid(size(tape(t)%patches)))
        else
          allocate(tape(t)%hlist(f)%varid(1))
        end if
      end if
      !
      ! If this field will be put out as columns then get column names for field
      !
      if (associated(tape(t)%patches)) then
        num_patches = size(tape(t)%patches)
        fldname = strip_suffix(tape(t)%hlist(f)%field%name)
        do i = 1, num_patches
          fname_tmp = trim(fldname)
          call tape(t)%patches(i)%field_name(fname_tmp)
          ierr = pio_inq_varid(tape(t)%File, trim(fname_tmp), tape(t)%hlist(f)%varid(i))
          call cam_pio_handle_error(ierr, 'H_INQUIRE: Error getting ID for '//trim(fname_tmp))
          ierr = pio_get_att(tape(t)%File, tape(t)%hlist(f)%varid(i), 'basename', basename)
          call cam_pio_handle_error(ierr, 'H_INQUIRE: Error getting basename for '//trim(fname_tmp))
          if (trim(fldname) /= trim(basename)) then
            call endrun('H_INQUIRE: basename ('//trim(basename)//') does not match fldname ('//trim(fldname)//')')
          end if
        end do
      else
        fldname = tape(t)%hlist(f)%field%name
        ierr = pio_inq_varid(tape(t)%File, trim(fldname), tape(t)%hlist(f)%varid(1))
        call cam_pio_handle_error(ierr, 'H_INQUIRE: Error getting ID for '//trim(fldname))
      end if
      if(tape(t)%hlist(f)%field%numlev>1) then
        ierr = pio_inq_attlen(tape(t)%File,tape(t)%hlist(f)%varid(1),'mdims', mdimsize)
        if(.not. associated(tape(t)%hlist(f)%field%mdims)) then
          allocate(tape(t)%hlist(f)%field%mdims(mdimsize))
        end if
        ierr=pio_get_att(tape(t)%File,tape(t)%hlist(f)%varid(1),'mdims', &
             tape(t)%hlist(f)%field%mdims(1:mdimsize))
        if(mdimsize>maxvarmdims) maxvarmdims=mdimsize
      end if

    end do

    if(masterproc) then
      write(iulog,*)'H_INQUIRE: Successfully opened netcdf file '
    end if

    return
  end subroutine h_inquire

  !#######################################################################

  subroutine add_default (name, tindex, flag)
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Add a field to the default "on" list for a given history file
    !
    ! Method:
    !
    !-----------------------------------------------------------------------
    !
    ! Arguments
    !
    character(len=*), intent(in) :: name  ! field name
    character(len=1), intent(in) :: flag  ! averaging flag

    integer, intent(in) :: tindex         ! history tape index
    !
    ! Local workspace
    !
    integer :: t            ! file index
    type(master_entry), pointer :: listentry

    if (htapes_defined) then
      call endrun ('ADD_DEFAULT: Attempt to add hist default '//trim(name)//' after history files set')
    end if
    !
    ! Check validity of input arguments
    !
    if (tindex > ptapes) then
      write(iulog,*)'ADD_DEFAULT: tape index=', tindex, ' is too big'
      call endrun
    end if

    ! Add to IC file if tindex = 0, reset to ptapes
    if (tindex == 0) then
      t = ptapes
      if ( .not. is_initfile(file_index=t) ) return
    else
      t = tindex
    end if

    if (verify(flag, HIST_AVG_FLAGS) /= 0) then
      call endrun ('ADD_DEFAULT: unknown averaging flag='//flag)
    end if
    !
    ! Look through master list for input field name.  When found, set active
    ! flag for that tape to true.  Also set averaging flag if told to use other
    ! than default.
    !
    listentry => get_entry_by_name(masterlinkedlist, trim(name))
    if(.not.associated(listentry)) then
      call endrun ('ADD_DEFAULT: field = "'//trim(name)//'" not found')
    end if
    listentry%actflag(t) = .true.
    if (flag /= ' ') then
      listentry%avgflag(t) = flag
      call AvgflagToString(flag, listentry%time_op(t))
    end if

    return
  end subroutine add_default

  !#######################################################################

  subroutine h_override (t)
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Override default history tape contents for a specific tape
    !
    ! Method: Copy the flag into the master field list
    !
    !-----------------------------------------------------------------------
    !
    ! Arguments
    !
    integer, intent(in) :: t         ! history tape index
    !
    ! Local workspace
    !
    character(len=1) :: avgflg       ! lcl equiv of avgflag_pertape(t) (to address xlf90 compiler bug)

    type(master_entry), pointer :: listentry


    avgflg = avgflag_pertape(t)


    listentry=>masterlinkedlist
    do while(associated(listentry))
      call AvgflagToString(avgflg, listentry%time_op(t))
      listentry%avgflag(t) = avgflag_pertape(t)
      listentry=>listentry%next_entry
    end do

  end subroutine h_override

  !#######################################################################

  subroutine h_define (t, restart)
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Define contents of history file t
    !
    ! Method: Issue the required netcdf wrapper calls to define the history file contents
    !
    !-----------------------------------------------------------------------
     use phys_control,    only: phys_getopts
    use cam_grid_support, only: cam_grid_header_info_t
    use cam_grid_support, only: cam_grid_write_attr, cam_grid_write_var
    use time_manager,     only: get_step_size, get_ref_date, timemgr_get_calendar_cf
    use cam_abortutils,   only: endrun
    use cam_pio_utils,    only: vdesc_ptr, cam_pio_handle_error, cam_pio_def_dim
    use cam_pio_utils,    only: cam_pio_createfile, cam_pio_def_var
    use sat_hist,         only: sat_hist_define

    !-----------------------------------------------------------------------

    !
    ! Input arguments
    !
    integer, intent(in) :: t   ! tape index
    logical, intent(in) :: restart
    !
    ! Local workspace
    !
    integer :: i, j            ! longitude, latitude indices
    integer :: grd             ! indices for looping through grids
    integer :: f               ! field index
    integer :: ncreal          ! real data type for output
    integer :: dtime           ! timestep size
    integer :: sec_nhtfrq      ! nhtfrq converted to seconds
    integer :: ndbase = 0      ! days component of base time
    integer :: nsbase = 0      ! seconds component of base time
    integer :: nbdate          ! base date in yyyymmdd format
    integer :: nbsec           ! time of day component of base date [seconds]
    integer :: yr, mon, day    ! year, month, day components of a date

    character(len=max_chars) :: str       ! character temporary
    character(len=max_chars) :: fname_tmp ! local copy of field name
    character(len=max_chars) :: calendar  ! Calendar type
    character(len=max_chars) :: cell_methods ! For cell_methods attribute
    character(len=16)        :: time_per_freq
    character(len=128)       :: errormsg

    integer :: ret                        ! function return value

    !
    ! netcdf dimensions
    !
    integer :: chardim            ! character dimension id
    integer :: dimenchar(2)       ! character dimension ids
    integer :: nacsdims(2)        ! dimension ids for nacs (used in restart file)
    integer :: bnddim             ! bounds dimension id
    integer :: timdim             ! unlimited dimension id

    integer :: dimindex(8)        ! dimension ids for variable declaration
    integer :: dimids_tmp(8)      ! dimension ids for variable declaration

    !
    ! netcdf variables
    !
    ! A structure to hold the horizontal dimension and coordinate info
    type(cam_grid_header_info_t), allocatable :: header_info(:)
    ! For satellite files and column output
    type(vdesc_ptr), allocatable :: latvar(:)    ! latitude variable ids
    type(vdesc_ptr), allocatable :: lonvar(:)    ! longitude variable ids

    type(var_desc_t), pointer        :: varid => NULL() ! temporary variable descriptor
    integer                          :: num_hdims, fdims
    integer                          :: num_patches ! How many entries for a field on this tape?
    integer,          pointer        :: mdims(:) => NULL()
    integer                          :: mdimsize
    integer                          :: ierr
    integer,          allocatable    :: mdimids(:)
    integer                          :: amode
    logical                          :: interpolate
    logical                          :: patch_output
    integer                          :: cam_snapshot_before_num
    integer                          :: cam_snapshot_after_num
    character(len=32)                :: cam_take_snapshot_before
    character(len=32)                :: cam_take_snapshot_after

    call phys_getopts(cam_take_snapshot_before_out= cam_take_snapshot_before, &
                      cam_take_snapshot_after_out = cam_take_snapshot_after,  &
                      cam_snapshot_before_num_out = cam_snapshot_before_num,  &
                      cam_snapshot_after_num_out  = cam_snapshot_after_num)

    if(restart) then
      tape => restarthistory_tape
      if(masterproc) write(iulog,*)'Opening netcdf history restart file ', trim(hrestpath(t))
    else
      tape => history_tape
      if(masterproc) write(iulog,*)'Opening netcdf history file ', trim(nhfil(t))
    end if

    amode = PIO_CLOBBER

    if(restart) then
      call cam_pio_createfile (tape(t)%File, hrestpath(t), amode)
    else
      call cam_pio_createfile (tape(t)%File, nhfil(t), amode)
    end if
    if(is_satfile(t)) then
      interpolate = .false. ! !!XXgoldyXX: Do we ever want to support this?
      patch_output = .false.
      call cam_pio_def_dim(tape(t)%File, 'ncol', pio_unlimited, timdim)
      call cam_pio_def_dim(tape(t)%File, 'nbnd', 2, bnddim)

      allocate(latvar(1), lonvar(1))
      allocate(latvar(1)%vd, lonvar(1)%vd)
      call cam_pio_def_var(tape(t)%File, 'lat', pio_double, (/timdim/),       &
           latvar(1)%vd)
      ierr=pio_put_att (tape(t)%File, latvar(1)%vd, 'long_name', 'latitude')
      ierr=pio_put_att (tape(t)%File, latvar(1)%vd, 'units', 'degrees_north')

      call cam_pio_def_var(tape(t)%File, 'lon', pio_double, (/timdim/),       &
           lonvar(1)%vd)
      ierr=pio_put_att (tape(t)%File, lonvar(1)%vd,'long_name','longitude')
      ierr=pio_put_att (tape(t)%File, lonvar(1)%vd,'units','degrees_east')

    else
      !
      ! Setup netcdf file - create the dimensions of lat,lon,time,level
      !
      ! interpolate is only supported for unstructured dycores
      interpolate = (interpolate_output(t) .and. (.not. restart))
      patch_output = (associated(tape(t)%patches) .and. (.not. restart))

      ! First define the horizontal grid dims
      ! Interpolation is special in that we ignore the native grids
      if(interpolate) then
        allocate(header_info(1))
        call cam_grid_write_attr(tape(t)%File, interpolate_info(t)%grid_id, header_info(1))
      else if (patch_output) then
        ! We are doing patch (column) output
        if (allocated(header_info)) then
          ! We shouldn't have any header_info yet
          call endrun('H_DEFINE: header_info should not be allocated for patch output')
        end if
        do i = 1, size(tape(t)%patches)
          call tape(t)%patches(i)%write_attrs(tape(t)%File)
        end do
      else
        allocate(header_info(size(tape(t)%grid_ids)))
        do i = 1, size(tape(t)%grid_ids)
          call cam_grid_write_attr(tape(t)%File, tape(t)%grid_ids(i), header_info(i))
        end do
      end if   ! interpolate

      ! Define the unlimited time dim
      call cam_pio_def_dim(tape(t)%File, 'time', pio_unlimited, timdim)
      call cam_pio_def_dim(tape(t)%File, 'nbnd', 2, bnddim, existOK=.true.)
      call cam_pio_def_dim(tape(t)%File, 'chars', 8, chardim)
    end if   ! is satfile

    ! Store snapshot location
    if (t == cam_snapshot_before_num) then
       ierr=pio_put_att(tape(t)%File, PIO_GLOBAL, 'cam_snapshot_before',      &
            trim(cam_take_snapshot_before))
    end if
    if (t == cam_snapshot_after_num) then
       ierr=pio_put_att(tape(t)%File, PIO_GLOBAL, 'cam_snapshot_after',       &
            trim(cam_take_snapshot_after))
    end if

    ! Populate the history coordinate (well, mdims anyway) attributes
    ! This routine also allocates the mdimids array
    call write_hist_coord_attrs(tape(t)%File, bnddim, mdimids, restart)

    call get_ref_date(yr, mon, day, nbsec)
    nbdate = yr*10000 + mon*100 + day
    ierr=pio_def_var (tape(t)%File,'time',pio_double,(/timdim/),tape(t)%timeid)
    ierr=pio_put_att (tape(t)%File, tape(t)%timeid, 'long_name', 'time')
    str = 'days since ' // date2yyyymmdd(nbdate) // ' ' // sec2hms(nbsec)
    ierr=pio_put_att (tape(t)%File, tape(t)%timeid, 'units', trim(str))

    calendar = timemgr_get_calendar_cf()
    ierr=pio_put_att (tape(t)%File, tape(t)%timeid, 'calendar', trim(calendar))


    ierr=pio_def_var (tape(t)%File,'date    ',pio_int,(/timdim/),tape(t)%dateid)
    str = 'current date (YYYYMMDD)'
    ierr=pio_put_att (tape(t)%File, tape(t)%dateid, 'long_name', trim(str))


    ierr=pio_def_var (tape(t)%File,'datesec ',pio_int,(/timdim/), tape(t)%datesecid)
    str = 'current seconds of current date'
    ierr=pio_put_att (tape(t)%File, tape(t)%datesecid, 'long_name', trim(str))

    !
    ! Character header information
    !
    str = 'CF-1.0'
    ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'Conventions', trim(str))
    ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'source', 'CAM')
#if ( defined BFB_CAM_SCAM_IOP )
    ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'CAM_GENERATED_FORCING','create SCAM IOP dataset')
#endif
    ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'case',caseid)
    ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'logname',logname)
    ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'host', host)

! Put these back in when they are filled properly
!    ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'title',ctitle)
!    ierr= pio_put_att (tape(t)%File, PIO_GLOBAL, 'Version', &
!         '$Name$')
!    ierr= pio_put_att (tape(t)%File, PIO_GLOBAL, 'revision_Id', &
!         '$Id$')

    ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'initial_file', ncdata)
    ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'topography_file', bnd_topo)
    if (len_trim(model_doi_url) > 0) then
       ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'model_doi_url', model_doi_url)
    end if

    ! Determine what time period frequency is being output for each file
    ! Note that nhtfrq is now in timesteps

    sec_nhtfrq = nhtfrq(t)

    ! If nhtfrq is in hours, convert to seconds
    if (nhtfrq(t) < 0) then
      sec_nhtfrq = abs(nhtfrq(t))*3600
    end if

    dtime = get_step_size()
    if (sec_nhtfrq == 0) then                                !month
      time_per_freq = 'month_1'
    else if (mod(sec_nhtfrq*dtime,86400) == 0) then          ! day
      write(time_per_freq,999) 'day_',sec_nhtfrq*dtime/86400
    else if (mod(sec_nhtfrq*dtime,3600) == 0) then           ! hour
      write(time_per_freq,999) 'hour_',(sec_nhtfrq*dtime)/3600
    else if (mod(sec_nhtfrq*dtime,60) == 0) then             ! minute
      write(time_per_freq,999) 'minute_',(sec_nhtfrq*dtime)/60
    else                                                     ! second
      write(time_per_freq,999) 'second_',sec_nhtfrq*dtime
    end if
999 format(a,i0)

    ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'time_period_freq', trim(time_per_freq))

    if(.not. is_satfile(t)) then

      ierr=pio_put_att (tape(t)%File, tape(t)%timeid, 'bounds', 'time_bnds')

      ierr=pio_def_var (tape(t)%File,'time_bnds',pio_double,(/bnddim,timdim/),tape(t)%tbndid)
      ierr=pio_put_att (tape(t)%File, tape(t)%tbndid, 'long_name', 'time interval endpoints')
      !
      ! Character
      !
      dimenchar(1) = chardim
      dimenchar(2) = timdim
      ierr=pio_def_var (tape(t)%File,'date_written',PIO_CHAR,dimenchar, tape(t)%date_writtenid)
      ierr=pio_def_var (tape(t)%File,'time_written',PIO_CHAR,dimenchar, tape(t)%time_writtenid)
      !
      ! Integer Header
      !

      ierr=pio_def_var (tape(t)%File,'ndbase',PIO_INT,tape(t)%ndbaseid)
      str = 'base day'
      ierr=pio_put_att (tape(t)%File, tape(t)%ndbaseid, 'long_name', trim(str))

      ierr=pio_def_var (tape(t)%File,'nsbase',PIO_INT,tape(t)%nsbaseid)
      str = 'seconds of base day'
      ierr=pio_put_att (tape(t)%File, tape(t)%nsbaseid, 'long_name', trim(str))

      ierr=pio_def_var (tape(t)%File,'nbdate',PIO_INT,tape(t)%nbdateid)
      str = 'base date (YYYYMMDD)'
      ierr=pio_put_att (tape(t)%File, tape(t)%nbdateid, 'long_name', trim(str))

#if ( defined BFB_CAM_SCAM_IOP )
      ierr=pio_def_var (tape(t)%File,'bdate',PIO_INT,tape(t)%bdateid)
      str = 'base date (YYYYMMDD)'
      ierr=pio_put_att (tape(t)%File, tape(t)%bdateid, 'long_name', trim(str))
#endif
      ierr=pio_def_var (tape(t)%File,'nbsec',PIO_INT,tape(t)%nbsecid)
      str = 'seconds of base date'
      ierr=pio_put_att (tape(t)%File, tape(t)%nbsecid, 'long_name', trim(str))

      ierr=pio_def_var (tape(t)%File,'mdt',PIO_INT,tape(t)%mdtid)
      ierr=pio_put_att (tape(t)%File, tape(t)%mdtid, 'long_name', 'timestep')
      ierr=pio_put_att (tape(t)%File, tape(t)%mdtid, 'units', 's')

      !
      ! Create variables for model timing and header information
      !

      ierr=pio_def_var (tape(t)%File,'ndcur   ',pio_int,(/timdim/),tape(t)%ndcurid)
      str = 'current day (from base day)'
      ierr=pio_put_att (tape(t)%File, tape(t)%ndcurid, 'long_name', trim(str))

      ierr=pio_def_var (tape(t)%File,'nscur   ',pio_int,(/timdim/),tape(t)%nscurid)
      str = 'current seconds of current day'
      ierr=pio_put_att (tape(t)%File, tape(t)%nscurid, 'long_name', trim(str))


      if (.not. is_initfile(file_index=t)) then
        ! Don't write the GHG/Solar forcing data to the IC file.
        ierr=pio_def_var (tape(t)%File,'co2vmr  ',pio_double,(/timdim/),tape(t)%co2vmrid)
        str = 'co2 volume mixing ratio'
        ierr=pio_put_att (tape(t)%File, tape(t)%co2vmrid, 'long_name', trim(str))

        ierr=pio_def_var (tape(t)%File,'ch4vmr  ',pio_double,(/timdim/),tape(t)%ch4vmrid)
        str = 'ch4 volume mixing ratio'
        ierr=pio_put_att (tape(t)%File, tape(t)%ch4vmrid, 'long_name', trim(str))

        ierr=pio_def_var (tape(t)%File,'n2ovmr  ',pio_double,(/timdim/),tape(t)%n2ovmrid)
        str = 'n2o volume mixing ratio'
        ierr=pio_put_att (tape(t)%File, tape(t)%n2ovmrid, 'long_name', trim(str))

        ierr=pio_def_var (tape(t)%File,'f11vmr  ',pio_double,(/timdim/),tape(t)%f11vmrid)
        str = 'f11 volume mixing ratio'
        ierr=pio_put_att (tape(t)%File, tape(t)%f11vmrid, 'long_name', trim(str))

        ierr=pio_def_var (tape(t)%File,'f12vmr  ',pio_double,(/timdim/),tape(t)%f12vmrid)
        str = 'f12 volume mixing ratio'
        ierr=pio_put_att (tape(t)%File, tape(t)%f12vmrid, 'long_name', trim(str))

        ierr=pio_def_var (tape(t)%File,'sol_tsi ',pio_double,(/timdim/),tape(t)%sol_tsiid)
        str = 'total solar irradiance'
        ierr=pio_put_att (tape(t)%File, tape(t)%sol_tsiid, 'long_name', trim(str))
        str = 'W/m2'
        ierr=pio_put_att (tape(t)%File, tape(t)%sol_tsiid, 'units', trim(str))

        if (solar_parms_on) then
          ! solar / geomagetic activity indices...
          ierr=pio_def_var (tape(t)%File,'f107',pio_double,(/timdim/),tape(t)%f107id)
          str = '10.7 cm solar radio flux (F10.7)'
          ierr=pio_put_att (tape(t)%File, tape(t)%f107id, 'long_name', trim(str))
          str = '10^-22 W m^-2 Hz^-1'
          ierr=pio_put_att (tape(t)%File, tape(t)%f107id, 'units', trim(str))

          ierr=pio_def_var (tape(t)%File,'f107a',pio_double,(/timdim/),tape(t)%f107aid)
          str = '81-day centered mean of 10.7 cm solar radio flux (F10.7)'
          ierr=pio_put_att (tape(t)%File, tape(t)%f107aid, 'long_name', trim(str))

          ierr=pio_def_var (tape(t)%File,'f107p',pio_double,(/timdim/),tape(t)%f107pid)
          str = 'Pervious day 10.7 cm solar radio flux (F10.7)'
          ierr=pio_put_att (tape(t)%File, tape(t)%f107pid, 'long_name', trim(str))

          ierr=pio_def_var (tape(t)%File,'kp',pio_double,(/timdim/),tape(t)%kpid)
          str = 'Daily planetary K geomagnetic index'
          ierr=pio_put_att (tape(t)%File, tape(t)%kpid, 'long_name', trim(str))

          ierr=pio_def_var (tape(t)%File,'ap',pio_double,(/timdim/),tape(t)%apid)
          str = 'Daily planetary A geomagnetic index'
          ierr=pio_put_att (tape(t)%File, tape(t)%apid, 'long_name', trim(str))
        endif
        if (solar_wind_on) then

          ierr=pio_def_var (tape(t)%File,'byimf',pio_double,(/timdim/),tape(t)%byimfid)
          str = 'Y component of the interplanetary magnetic field'
          ierr=pio_put_att (tape(t)%File, tape(t)%byimfid, 'long_name', trim(str))
          str = 'nT'
          ierr=pio_put_att (tape(t)%File, tape(t)%byimfid, 'units', trim(str))

          ierr=pio_def_var (tape(t)%File,'bzimf',pio_double,(/timdim/),tape(t)%bzimfid)
          str = 'Z component of the interplanetary magnetic field'
          ierr=pio_put_att (tape(t)%File, tape(t)%bzimfid, 'long_name', trim(str))
          str = 'nT'
          ierr=pio_put_att (tape(t)%File, tape(t)%bzimfid, 'units', trim(str))

          ierr=pio_def_var (tape(t)%File,'swvel',pio_double,(/timdim/),tape(t)%swvelid)
          str = 'Solar wind speed'
          ierr=pio_put_att (tape(t)%File, tape(t)%swvelid, 'long_name', trim(str))
          str = 'km/sec'
          ierr=pio_put_att (tape(t)%File, tape(t)%swvelid, 'units', trim(str))

          ierr=pio_def_var (tape(t)%File,'swden',pio_double,(/timdim/),tape(t)%swdenid)
          str = 'Solar wind ion number density'
          ierr=pio_put_att (tape(t)%File, tape(t)%swdenid, 'long_name', trim(str))
          str = 'cm-3'
          ierr=pio_put_att (tape(t)%File, tape(t)%swdenid, 'units', trim(str))

        endif
      end if


#if ( defined BFB_CAM_SCAM_IOP )
      ierr=pio_def_var (tape(t)%File,'tsec ',pio_int,(/timdim/), tape(t)%tsecid)
      str = 'current seconds of current date needed for scam'
      ierr=pio_put_att (tape(t)%File, tape(t)%tsecid, 'long_name', trim(str))
#endif
      ierr=pio_def_var (tape(t)%File,'nsteph  ',pio_int,(/timdim/),tape(t)%nstephid)
      str = 'current timestep'
      ierr=pio_put_att (tape(t)%File, tape(t)%nstephid, 'long_name', trim(str))
    end if ! .not. is_satfile

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !
    ! Create variables and attributes for field list
    !
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    do f = 1, nflds(t)

      !! Collect some field properties
      call AvgflagToString(tape(t)%hlist(f)%avgflag, tape(t)%hlist(f)%time_op)

      if ((tape(t)%hlist(f)%hwrt_prec == 8) .or. restart) then
        ncreal = pio_double
      else
        ncreal = pio_real
      end if

      if(associated(tape(t)%hlist(f)%field%mdims)) then
        mdims => tape(t)%hlist(f)%field%mdims
        mdimsize = size(mdims)
      else if(tape(t)%hlist(f)%field%numlev > 1) then
        call endrun('mdims not defined for variable '//trim(tape(t)%hlist(f)%field%name))
      else
        mdimsize=0
      end if

      ! num_patches will loop through the number of patches (or just one
      !             for the whole grid) for this field for this tape
      if (patch_output) then
        num_patches = size(tape(t)%patches)
      else
        num_patches = 1
      end if
      if(.not.associated(tape(t)%hlist(f)%varid)) then
        allocate(tape(t)%hlist(f)%varid(num_patches))
      end if
      fname_tmp = strip_suffix(tape(t)%hlist(f)%field%name)

      if(is_satfile(t)) then
        num_hdims=0
        nfils(t)=1
        call sat_hist_define(tape(t)%File)
      else if (interpolate) then
        ! Interpolate can't use normal grid code since we are forcing fields
        ! to use interpolate decomp
        if (.not. allocated(header_info)) then
          ! Safety check
          call endrun('h_define: header_info not allocated')
        end if
        num_hdims = 2
        do i = 1, num_hdims
          dimindex(i) = header_info(1)%get_hdimid(i)
          nacsdims(i) = header_info(1)%get_hdimid(i)
        end do
      else if (patch_output) then
        ! All patches for this variable should be on the same grid
        num_hdims = tape(t)%patches(1)%num_hdims(tape(t)%hlist(f)%field%decomp_type)
      else
        ! Normal grid output
        ! Find appropriate grid in header_info
        if (.not. allocated(header_info)) then
          ! Safety check
          call endrun('h_define: header_info not allocated')
        end if
        grd = -1
        do i = 1, size(header_info)
          if (header_info(i)%get_gridid() == tape(t)%hlist(f)%field%decomp_type) then
            grd = i
            exit
          end if
        end do
        if (grd < 0) then
          write(errormsg, '(a,i0,2a)') 'grid, ',tape(t)%hlist(f)%field%decomp_type,', not found for ',trim(fname_tmp)
          call endrun('H_DEFINE: '//errormsg)
        end if
        num_hdims = header_info(grd)%num_hdims()
        do i = 1, num_hdims
          dimindex(i) = header_info(grd)%get_hdimid(i)
          nacsdims(i) = header_info(grd)%get_hdimid(i)
        end do
      end if     ! is_satfile

      !
      !  Create variables and atributes for fields written out as columns
      !

      do i = 1, num_patches
        fname_tmp = strip_suffix(tape(t)%hlist(f)%field%name)
        varid => tape(t)%hlist(f)%varid(i)
        dimids_tmp = dimindex
        ! Figure the dimension ID array for this field
        ! We have defined the horizontal grid dimensions in dimindex
        fdims = num_hdims
        do j = 1, mdimsize
          fdims = fdims + 1
          dimids_tmp(fdims) = mdimids(mdims(j))
        end do
        if(.not. restart) then
          ! Only add time dimension if this is not a restart history tape
          fdims = fdims + 1
          dimids_tmp(fdims) = timdim
        end if
        if (patch_output) then
          ! For patch output, we need new dimension IDs and a different name
          call tape(t)%patches(i)%get_var_data(fname_tmp,                     &
               dimids_tmp(1:fdims), tape(t)%hlist(f)%field%decomp_type)
        end if
        ! Define the variable
        call cam_pio_def_var(tape(t)%File, trim(fname_tmp), ncreal,           &
             dimids_tmp(1:fdims), varid)
        if (mdimsize > 0) then
          ierr = pio_put_att(tape(t)%File, varid, 'mdims', mdims(1:mdimsize))
          call cam_pio_handle_error(ierr, 'h_define: cannot define mdims for '//trim(fname_tmp))
        end if
        str = tape(t)%hlist(f)%field%sampling_seq
        if (len_trim(str) > 0) then
          ierr = pio_put_att(tape(t)%File, varid, 'Sampling_Sequence', trim(str))
          call cam_pio_handle_error(ierr, 'h_define: cannot define Sampling_Sequence for '//trim(fname_tmp))
        end if

        if (tape(t)%hlist(f)%field%flag_xyfill) then
          ! Add both _FillValue and missing_value to cover expectations
          !     of various applications.
          ! The attribute type must match the data type.
          if ((tape(t)%hlist(f)%hwrt_prec == 8) .or. restart) then
            ierr = pio_put_att(tape(t)%File, varid, '_FillValue',             &
                 tape(t)%hlist(f)%field%fillvalue)
            call cam_pio_handle_error(ierr,                                   &
                 'h_define: cannot define _FillValue for '//trim(fname_tmp))
            ierr = pio_put_att(tape(t)%File, varid, 'missing_value',          &
                 tape(t)%hlist(f)%field%fillvalue)
            call cam_pio_handle_error(ierr,                                   &
                 'h_define: cannot define missing_value for '//trim(fname_tmp))
          else
            ierr = pio_put_att(tape(t)%File, varid, '_FillValue',             &
                 REAL(tape(t)%hlist(f)%field%fillvalue,r4))
            call cam_pio_handle_error(ierr,                                   &
                 'h_define: cannot define _FillValue for '//trim(fname_tmp))
            ierr = pio_put_att(tape(t)%File, varid, 'missing_value',          &
                 REAL(tape(t)%hlist(f)%field%fillvalue,r4))
            call cam_pio_handle_error(ierr,                                   &
                 'h_define: cannot define missing_value for '//trim(fname_tmp))
          end if
        end if

        str = tape(t)%hlist(f)%field%units
        if (len_trim(str) > 0) then
          ierr=pio_put_att (tape(t)%File, varid, 'units', trim(str))
          call cam_pio_handle_error(ierr,                                     &
               'h_define: cannot define units for '//trim(fname_tmp))
        end if

        str = tape(t)%hlist(f)%field%mixing_ratio
        if (len_trim(str) > 0) then
          ierr=pio_put_att (tape(t)%File, varid, 'mixing_ratio', trim(str))
          call cam_pio_handle_error(ierr,                                     &
               'h_define: cannot define mixing_ratio for '//trim(fname_tmp))
        end if

        str = tape(t)%hlist(f)%field%long_name
        ierr=pio_put_att (tape(t)%File, varid, 'long_name', trim(str))
        call cam_pio_handle_error(ierr,                                       &
             'h_define: cannot define long_name for '//trim(fname_tmp))

        ! Assign field attributes defining valid levels and averaging info

        cell_methods = ''
        if (len_trim(tape(t)%hlist(f)%field%cell_methods) > 0) then
          if (len_trim(cell_methods) > 0) then
            cell_methods = trim(cell_methods)//' '//trim(tape(t)%hlist(f)%field%cell_methods)
          else
            cell_methods = trim(cell_methods)//trim(tape(t)%hlist(f)%field%cell_methods)
          end if
        end if
        ! Time cell methods is after field method because time averaging is
        ! applied later (just before output) than field method which is applied
        ! before outfld call.
        str = tape(t)%hlist(f)%time_op
        select case (str)
        case ('mean', 'maximum', 'minimum', 'standard_deviation')
          if (len_trim(cell_methods) > 0) then
            cell_methods = trim(cell_methods)//' '//'time: '//str
          else
            cell_methods = trim(cell_methods)//'time: '//str
          end if
        end select
        if (len_trim(cell_methods) > 0) then
          ierr = pio_put_att(tape(t)%File, varid, 'cell_methods', trim(cell_methods))
          call cam_pio_handle_error(ierr,                                     &
               'h_define: cannot define cell_methods for '//trim(fname_tmp))
        end if
        if (patch_output) then
          ierr = pio_put_att(tape(t)%File, varid, 'basename',                 &
               tape(t)%hlist(f)%field%name)
          call cam_pio_handle_error(ierr,                                     &
               'h_define: cannot define basename for '//trim(fname_tmp))
        end if

        if (restart) then
          ! For restart history files, we need to save accumulation counts
          fname_tmp = trim(fname_tmp)//'_nacs'
          if (.not. associated(tape(t)%hlist(f)%nacs_varid)) then
            allocate(tape(t)%hlist(f)%nacs_varid)
          end if
          if (size(tape(t)%hlist(f)%nacs, 1) > 1) then
            call cam_pio_def_var(tape(t)%File, trim(fname_tmp), pio_int,      &
                 nacsdims(1:num_hdims), tape(t)%hlist(f)%nacs_varid)
          else
            ! Save just one value representing all chunks
            call cam_pio_def_var(tape(t)%File, trim(fname_tmp), pio_int,      &
                 tape(t)%hlist(f)%nacs_varid)
          end if
          ! for standard deviation
          if (associated(tape(t)%hlist(f)%sbuf)) then
             fname_tmp = strip_suffix(tape(t)%hlist(f)%field%name)
             fname_tmp = trim(fname_tmp)//'_var'
             if ( .not.associated(tape(t)%hlist(f)%sbuf_varid)) then
                allocate(tape(t)%hlist(f)%sbuf_varid)
             endif
             call cam_pio_def_var(tape(t)%File, trim(fname_tmp), pio_double,      &
                  dimids_tmp(1:fdims), tape(t)%hlist(f)%sbuf_varid)
          endif
        end if
      end do ! Loop over output patches
    end do   ! Loop over fields
    !
    deallocate(mdimids)
    ret = pio_enddef(tape(t)%File)

    if(masterproc) then
      write(iulog,*)'H_DEFINE: Successfully opened netcdf file '
    endif
    !
    ! Write time-invariant portion of history header
    !
    if(.not. is_satfile(t)) then
      if(interpolate) then
        call cam_grid_write_var(tape(t)%File, interpolate_info(t)%grid_id)
      else if((.not. patch_output) .or. restart) then
        do i = 1, size(tape(t)%grid_ids)
          call cam_grid_write_var(tape(t)%File, tape(t)%grid_ids(i))
        end do
      else
        ! Patch output
        do i = 1, size(tape(t)%patches)
          call tape(t)%patches(i)%write_vals(tape(t)%File)
        end do
      end if ! interpolate
      if (allocated(lonvar)) then
        deallocate(lonvar)
      end if
      if (allocated(latvar)) then
        deallocate(latvar)
      end if

      dtime = get_step_size()
      ierr = pio_put_var(tape(t)%File, tape(t)%mdtid, (/dtime/))
      call cam_pio_handle_error(ierr, 'h_define: cannot put mdt')
      !
      ! Model date info
      !
      ierr = pio_put_var(tape(t)%File, tape(t)%ndbaseid, (/ndbase/))
      call cam_pio_handle_error(ierr, 'h_define: cannot put ndbase')
      ierr = pio_put_var(tape(t)%File, tape(t)%nsbaseid, (/nsbase/))
      call cam_pio_handle_error(ierr, 'h_define: cannot put nsbase')

      ierr = pio_put_var(tape(t)%File, tape(t)%nbdateid, (/nbdate/))
      call cam_pio_handle_error(ierr, 'h_define: cannot put nbdate')
#if ( defined BFB_CAM_SCAM_IOP )
      ierr = pio_put_var(tape(t)%File, tape(t)%bdateid, (/nbdate/))
      call cam_pio_handle_error(ierr, 'h_define: cannot put bdate')
#endif
      ierr = pio_put_var(tape(t)%File, tape(t)%nbsecid, (/nbsec/))
      call cam_pio_handle_error(ierr, 'h_define: cannot put nbsec')
      !
      ! Reduced grid info
      !

    end if ! .not. is_satfile

    if (allocated(header_info)) then
      do i = 1, size(header_info)
        call header_info(i)%deallocate()
      end do
      deallocate(header_info)
    end if

    ! Write the mdim variable data
    call write_hist_coord_vars(tape(t)%File, restart)

  end subroutine h_define

  !#######################################################################

  subroutine h_normalize (f, t)

    use cam_history_support, only: dim_index_2d

    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Normalize fields on a history file by the number of accumulations
    !
    ! Method: Loop over fields on the tape.  Need averaging flag and number of
    !         accumulations to perform normalization.
    !
    !-----------------------------------------------------------------------
    !
    ! Input arguments
    !
    integer, intent(in) :: f       ! field index
    integer, intent(in) :: t       ! tape index
    !
    ! Local workspace
    !
    type (dim_index_2d) :: dimind  ! 2-D dimension index
    integer     :: c               ! chunk (or lat) index
    integer     :: ib, ie    ! beginning and ending indices of first dimension
    integer     :: jb, je    ! beginning and ending indices of second dimension
    integer     :: begdim3, enddim3 ! Chunk or block bounds
    integer     :: k         ! level
    integer     :: i, ii
    real(r8) :: variance, tmpfill

    logical     :: flag_xyfill ! non-applicable xy points flagged with fillvalue
    character*1 :: avgflag     ! averaging flag

    call t_startf ('h_normalize')

    call tape(t)%hlist(f)%field%get_bounds(3, begdim3, enddim3)

    !
    ! normalize by number of accumulations for averaged case
    !
    flag_xyfill = tape(t)%hlist(f)%field%flag_xyfill
    avgflag = tape(t)%hlist(f)%avgflag

    do c = begdim3, enddim3
      dimind = tape(t)%hlist(f)%field%get_dims(c)

      ib = dimind%beg1
      ie = dimind%end1
      jb = dimind%beg2
      je = dimind%end2

      if (flag_xyfill) then
        do k = jb, je
          where (tape(t)%hlist(f)%nacs(ib:ie, c) == 0)
            tape(t)%hlist(f)%hbuf(ib:ie,k,c) = tape(t)%hlist(f)%field%fillvalue
          endwhere
        end do
      end if

      if (avgflag == 'A' .or. avgflag == 'B' .or. avgflag == 'L') then
        if (size(tape(t)%hlist(f)%nacs, 1) > 1) then
          do k = jb, je
            where (tape(t)%hlist(f)%nacs(ib:ie,c) /= 0)
              tape(t)%hlist(f)%hbuf(ib:ie,k,c) = &
                   tape(t)%hlist(f)%hbuf(ib:ie,k,c) &
                   / tape(t)%hlist(f)%nacs(ib:ie,c)
            endwhere
          end do
        else if(tape(t)%hlist(f)%nacs(1,c) > 0) then
          do k=jb,je
            tape(t)%hlist(f)%hbuf(ib:ie,k,c) = &
                 tape(t)%hlist(f)%hbuf(ib:ie,k,c) &
                 / tape(t)%hlist(f)%nacs(1,c)
          end do
        end if
      end if
      if (avgflag == 'S') then
         ! standard deviation ...
         ! from http://www.johndcook.com/blog/standard_deviation/
         tmpfill = merge(tape(t)%hlist(f)%field%fillvalue,0._r8,flag_xyfill)
         do k=jb,je
            do i = ib,ie
               ii = merge(i,1,flag_xyfill)
               if (tape(t)%hlist(f)%nacs(ii,c) > 1) then
                  variance = tape(t)%hlist(f)%sbuf(i,k,c)/(tape(t)%hlist(f)%nacs(ii,c)-1)
                  tape(t)%hlist(f)%hbuf(i,k,c) = sqrt(variance)
               else
                  tape(t)%hlist(f)%hbuf(i,k,c) = tmpfill
               endif
            end do
         end do
      endif
    end do

    call t_stopf ('h_normalize')

    return
  end subroutine h_normalize

  !#######################################################################

  subroutine h_zero (f, t)
    use cam_history_support, only: dim_index_2d
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Zero out accumulation buffers for a tape
    !
    ! Method: Loop through fields on the tape
    !
    !-----------------------------------------------------------------------
    !
    integer, intent(in) :: f     ! field index
    integer, intent(in) :: t     ! tape index
    !
    ! Local workspace
    !
    type (dim_index_2d) :: dimind   ! 2-D dimension index
    integer :: c                    ! chunk index
    integer :: begdim3              ! on-node chunk or lat start index
    integer :: enddim3              ! on-node chunk or lat end index

    call t_startf ('h_zero')

    call tape(t)%hlist(f)%field%get_bounds(3, begdim3, enddim3)

    do c = begdim3, enddim3
      dimind = tape(t)%hlist(f)%field%get_dims(c)
      tape(t)%hlist(f)%hbuf(dimind%beg1:dimind%end1,dimind%beg2:dimind%end2,c)=0._r8
      if (associated(tape(t)%hlist(f)%sbuf)) then ! zero out variance buffer for standard deviation
         tape(t)%hlist(f)%sbuf(dimind%beg1:dimind%end1,dimind%beg2:dimind%end2,c)=0._r8
      endif
    end do
    tape(t)%hlist(f)%nacs(:,:) = 0

    call t_stopf ('h_zero')

    return
  end subroutine h_zero

  !#######################################################################

  subroutine dump_field (f, t, restart)
    use cam_history_support, only: history_patch_t, dim_index_3d
    use cam_grid_support,    only: cam_grid_write_dist_array, cam_grid_dimensions
    use cam_interp_mod,      only : write_interpolated

    ! Dummy arguments
    integer,     intent(in)    :: f
    integer,     intent(in)    :: t
    logical,     intent(in)    :: restart
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Write a variable to a history tape using PIO
    !          For restart tapes, also write the accumulation buffer (nacs)
    !
    !-----------------------------------------------------------------------
    ! Local variables
    integer                          :: ierr
    type(var_desc_t),      pointer   :: varid      ! PIO ID for var
    type(var_desc_t),      pointer   :: compid     ! PIO ID for vector comp.
    integer                          :: compind    ! index of vector comp.
    integer                          :: fdims(8)   ! Field file dim sizes
    integer                          :: frank      ! Field file rank
    integer                          :: nacsrank   ! Field file rank for nacs
    type(dim_index_3d)               :: dimind     ! 3-D dimension index
    integer                          :: adims(3)   ! Field array dim sizes
    integer                          :: nadims     ! # of used adims
    integer                          :: fdecomp
    integer                          :: num_patches
    integer                          :: mdimsize   ! Total # on-node elements
    logical                          :: interpolate
    logical                          :: patch_output
    type(history_patch_t), pointer   :: patchptr
    integer :: i

    interpolate = (interpolate_output(t) .and. (.not. restart))
    patch_output = (associated(tape(t)%patches) .and. (.not. restart))

    !!! Get the field's shape and decomposition

    ! Shape on disk
    call tape(t)%hlist(f)%field%get_shape(fdims, frank)

    ! Shape of array
    dimind = tape(t)%hlist(f)%field%get_dims()
    call dimind%dim_sizes(adims)
    if (adims(2) <= 1) then
      adims(2) = adims(3)
      nadims = 2
    else
      nadims = 3
    end if
    fdecomp = tape(t)%hlist(f)%field%decomp_type

    ! num_patches will loop through the number of patches (or just one
    !             for the whole grid) for this field for this tape
    if (patch_output) then
      num_patches = size(tape(t)%patches)
    else
      num_patches = 1
    end if

    do i = 1, num_patches
      varid => tape(t)%hlist(f)%varid(i)

      if (restart) then
        call pio_setframe(tape(t)%File, varid, int(-1,kind=PIO_OFFSET_KIND))
      else
        call pio_setframe(tape(t)%File, varid, int(max(1,nfils(t)),kind=PIO_OFFSET_KIND))
      end if
      if (patch_output) then
        ! We are outputting patches
        patchptr => tape(t)%patches(i)
        if (interpolate) then
          call endrun('dump_field: interpolate incompatible with regional output')
        end if
        call patchptr%write_var(tape(t)%File, fdecomp, adims(1:nadims),       &
             pio_double, tape(t)%hlist(f)%hbuf, varid)
      else
        ! We are doing output via the field's grid
        if (interpolate) then
          mdimsize = tape(t)%hlist(f)%field%enddim2 - tape(t)%hlist(f)%field%begdim2 + 1
          if (mdimsize == 0) then
            mdimsize = tape(t)%hlist(f)%field%numlev
          end if
          if (tape(t)%hlist(f)%field%meridional_complement > 0) then
            compind = tape(t)%hlist(f)%field%meridional_complement
            compid => tape(t)%hlist(compind)%varid(i)
            ! We didn't call set frame on the meridional complement field
            call pio_setframe(tape(t)%File, compid, int(max(1,nfils(t)),kind=PIO_OFFSET_KIND))
            call write_interpolated(tape(t)%File, varid, compid,              &
                 tape(t)%hlist(f)%hbuf, tape(t)%hlist(compind)%hbuf,          &
                 mdimsize, PIO_DOUBLE, fdecomp)
          else if (tape(t)%hlist(f)%field%zonal_complement > 0) then
            ! We don't want to double write so do nothing here
!            compind = tape(t)%hlist(f)%field%zonal_complement
!            compid => tape(t)%hlist(compind)%varid(i)
!            call write_interpolated(tape(t)%File, compid, varid,              &
!                 tape(t)%hlist(compind)%hbuf, tape(t)%hlist(f)%hbuf,          &
!                 mdimsize, PIO_DOUBLE, fdecomp)
          else
            ! Scalar field
            call write_interpolated(tape(t)%File, varid,                      &
                 tape(t)%hlist(f)%hbuf, mdimsize, PIO_DOUBLE, fdecomp)
          end if
        else if (nadims == 2) then
          ! Special case for 2D field (no levels) due to hbuf structure
          call cam_grid_write_dist_array(tape(t)%File, fdecomp,               &
               adims(1:nadims), fdims(1:frank), tape(t)%hlist(f)%hbuf(:,1,:), varid)
        else
          call cam_grid_write_dist_array(tape(t)%File, fdecomp, adims,        &
               fdims(1:frank), tape(t)%hlist(f)%hbuf, varid)
        end if
      end if
    end do
    !! write accumulation counter and variance to hist restart file
    if(restart) then
       if (associated(tape(t)%hlist(f)%sbuf) ) then
           ! write variance data to restart file for standard deviation calc
          if (nadims == 2) then
           ! Special case for 2D field (no levels) due to sbuf structure
             call cam_grid_write_dist_array(tape(t)%File, fdecomp, adims(1:nadims), &
                  fdims(1:frank), tape(t)%hlist(f)%sbuf(:,1,:), tape(t)%hlist(f)%sbuf_varid)
          else
             call cam_grid_write_dist_array(tape(t)%File, fdecomp, adims,        &
                  fdims(1:frank), tape(t)%hlist(f)%sbuf, tape(t)%hlist(f)%sbuf_varid)
          endif
       endif
     !! NACS
       if (size(tape(t)%hlist(f)%nacs, 1) > 1) then
          if (nadims > 2) then
             adims(2) = adims(3)
             nadims = 2
          end if
          call cam_grid_dimensions(fdecomp, fdims(1:2), nacsrank)
          call cam_grid_write_dist_array(tape(t)%File, fdecomp, adims(1:nadims), &
               fdims(1:nacsrank), tape(t)%hlist(f)%nacs, tape(t)%hlist(f)%nacs_varid)
       else
          ierr = pio_put_var(tape(t)%File, tape(t)%hlist(f)%nacs_varid,     &
               tape(t)%hlist(f)%nacs(:, tape(t)%hlist(f)%field%begdim3:tape(t)%hlist(f)%field%enddim3))
       end if
    end if

    return
  end subroutine dump_field

  !#######################################################################

  logical function write_inithist ()
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Set flags that will initiate dump to IC file when OUTFLD and
    ! WSHIST are called
    !
    !-----------------------------------------------------------------------
    !
    use time_manager, only: get_nstep, get_curr_date, get_step_size, is_last_step
    !
    ! Local workspace
    !
    integer :: yr, mon, day      ! year, month, and day components of
    ! a date
    integer :: nstep             ! current timestep number
    integer :: ncsec             ! current time of day [seconds]
    integer :: dtime             ! timestep size

    !-----------------------------------------------------------------------

    write_inithist  = .false.

    if(is_initfile()) then

      nstep = get_nstep()
      call get_curr_date(yr, mon, day, ncsec)

      if    (inithist == '6-HOURLY') then
        dtime  = get_step_size()
        write_inithist = nstep /= 0 .and. mod( nstep, nint((6._r8*3600._r8)/dtime) ) == 0
      elseif(inithist == 'DAILY'   ) then
        write_inithist = nstep /= 0 .and. ncsec == 0
      elseif(inithist == 'MONTHLY' ) then
        write_inithist = nstep /= 0 .and. ncsec == 0 .and. day == 1
      elseif(inithist == 'YEARLY'  ) then
        write_inithist = nstep /= 0 .and. ncsec == 0 .and. day == 1 .and. mon == 1
      elseif(inithist == 'CAMIOP'  ) then
        write_inithist = nstep == 0
      elseif(inithist == 'ENDOFRUN'  ) then
        write_inithist = nstep /= 0 .and. is_last_step()
      end if
    end if

    return
  end function write_inithist

  !#######################################################################

  subroutine wshist (rgnht_in)
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Driver routine to write fields on history tape t
    !
    !
    !-----------------------------------------------------------------------
    use time_manager,      only: get_nstep, get_curr_date, get_curr_time, get_step_size
    use chem_surfvals,     only: chem_surfvals_get, chem_surfvals_co2_rad
    use solar_irrad_data,  only: sol_tsi
    use sat_hist,          only: sat_hist_write
    use cam_interp_mod,    only: set_interp_hfile
    use cam_datetime_mod,  only: datetime
    use cam_pio_utils,     only: cam_pio_closefile

    logical, intent(in), optional :: rgnht_in(ptapes)
    !
    ! Local workspace
    !
    character(len=8) :: cdate  ! system date
    character(len=8) :: ctime  ! system time

    logical  :: rgnht(ptapes), restart
    integer t, f               ! tape, field indices
    integer start              ! starting index required by nf_put_vara
    integer count1             ! count values required by nf_put_vara
    integer startc(2)          ! start values required by nf_put_vara (character)
    integer countc(2)          ! count values required by nf_put_vara (character)
#ifdef HDEBUG
    !      integer begdim3
    !      integer enddim3
#endif

    integer :: yr, mon, day      ! year, month, and day components of a date
    integer :: nstep             ! current timestep number
    integer :: ncdate            ! current date in integer format [yyyymmdd]
    integer :: ncsec             ! current time of day [seconds]
    integer :: ndcur             ! day component of current time
    integer :: nscur             ! seconds component of current time
    real(r8) :: time             ! current time
    real(r8) :: tdata(2)         ! time interval boundaries
    character(len=max_string_len) :: fname ! Filename
    logical :: prev              ! Label file with previous date rather than current
    integer :: ierr
#if ( defined BFB_CAM_SCAM_IOP )
    integer :: tsec             ! day component of current time
    integer :: dtime            ! seconds component of current time
#endif

    if(present(rgnht_in)) then
      rgnht=rgnht_in
      restart=.true.
      tape => restarthistory_tape
    else
      rgnht=.false.
      restart=.false.
      tape => history_tape
    end if

    nstep = get_nstep()
    call get_curr_date(yr, mon, day, ncsec)
    ncdate = yr*10000 + mon*100 + day
    call get_curr_time(ndcur, nscur)
    !
    ! Write time-varying portion of history file header
    !
    do t=1,ptapes
      if (nflds(t) == 0 .or. (restart .and.(.not.rgnht(t)))) cycle
      !
      ! Check if this is the IC file and if it's time to write.
      ! Else, use "nhtfrq" to determine if it's time to write
      ! the other history files.
      !
      if((.not. restart) .or. rgnht(t)) then
        if( is_initfile(file_index=t) ) then
          hstwr(t) =  write_inithist()
          prev     = .false.
        else
          if (nhtfrq(t) == 0) then
            hstwr(t) = nstep /= 0 .and. day == 1 .and. ncsec == 0
            prev     = .true.
          else
            hstwr(t) = mod(nstep,nhtfrq(t)) == 0
            prev     = .false.
          end if
        end if
      end if
      if (hstwr(t) .or. (restart .and. rgnht(t))) then
        if(masterproc) then
          if(is_initfile(file_index=t)) then
            write(iulog,100) yr,mon,day,ncsec
100         format('WSHIST: writing time sample to Initial Conditions h-file', &
                 ' DATE=',i4.4,'/',i2.2,'/',i2.2,' NCSEC=',i6)
          else if(is_satfile(t)) then
            write(iulog,150) nfils(t),t,yr,mon,day,ncsec
150         format('WSHIST: writing sat columns ',i6,' to h-file ', &
                 i1,' DATE=',i4.4,'/',i2.2,'/',i2.2,' NCSEC=',i6)
          else if(hstwr(t)) then
            write(iulog,200) nfils(t),t,yr,mon,day,ncsec
200         format('WSHIST: writing time sample ',i3,' to h-file ', &
                 i1,' DATE=',i4.4,'/',i2.2,'/',i2.2,' NCSEC=',i6)
          else if(restart .and. rgnht(t)) then
            write(iulog,300) nfils(t),t,yr,mon,day,ncsec
300         format('WSHIST: writing history restart ',i3,' to hr-file ', &
                 i1,' DATE=',i4.4,'/',i2.2,'/',i2.2,' NCSEC=',i6)
          end if
          write(iulog,*)
        end if
        !
        ! Starting a new volume => define the metadata
        !
        if (nfils(t)==0 .or. (restart.and.rgnht(t))) then
          if(restart) then
            rhfilename_spec = '%c.cam' // trim(inst_suffix) // '.rh%t.%y-%m-%d-%s.nc'
            fname = interpret_filename_spec( rhfilename_spec, number=(t-1))
            hrestpath(t)=fname
          else if(is_initfile(file_index=t)) then
            fname = interpret_filename_spec( hfilename_spec(t) )
          else
            fname = interpret_filename_spec( hfilename_spec(t), number=(t-1), &
                 prev=prev )
          end if
          !
          ! Check that this new filename isn't the same as a previous or current filename
          !
          do f = 1, ptapes
            if (masterproc.and. trim(fname) == trim(nhfil(f)) )then
              write(iulog,*)'WSHIST: New filename same as old file = ', trim(fname)
              write(iulog,*)'Is there an error in your filename specifiers?'
              write(iulog,*)'hfilename_spec(', t, ') = ', hfilename_spec(t)
              if ( t /= f )then
                write(iulog,*)'hfilename_spec(', f, ') = ', hfilename_spec(f)
              end if
              call endrun
            end if
          end do
          if(.not. restart) then
            nhfil(t) = fname
            if(masterproc) write(iulog,*)'WSHIST: nhfil(',t,')=',trim(nhfil(t))
            cpath(t) = nhfil(t)
            if ( len_trim(nfpath(t)) == 0 ) nfpath(t) = cpath(t)
          end if
          call h_define (t, restart)
        end if

        if(is_satfile(t)) then
          call sat_hist_write( tape(t), nflds(t), nfils(t))
        else
          if(restart) then
            start=1
          else
            nfils(t) = nfils(t) + 1
            start = nfils(t)
          end if
          count1 = 1
          ! Setup interpolation data if history file is interpolated
          if (interpolate_output(t) .and. (.not. restart)) then
            call set_interp_hfile(t, interpolate_info)
          end if

          ierr = pio_put_var (tape(t)%File, tape(t)%ndcurid,(/start/), (/count1/),(/ndcur/))
          ierr = pio_put_var (tape(t)%File, tape(t)%nscurid,(/start/), (/count1/),(/nscur/))
          ierr = pio_put_var (tape(t)%File, tape(t)%dateid,(/start/), (/count1/),(/ncdate/))

          if (.not. is_initfile(file_index=t)) then
            ! Don't write the GHG/Solar forcing data to the IC file.
            ierr=pio_put_var (tape(t)%File, tape(t)%co2vmrid,(/start/), (/count1/),(/chem_surfvals_co2_rad(vmr_in=.true.)/))
            ierr=pio_put_var (tape(t)%File, tape(t)%ch4vmrid,(/start/), (/count1/),(/chem_surfvals_get('CH4VMR')/))
            ierr=pio_put_var (tape(t)%File, tape(t)%n2ovmrid,(/start/), (/count1/),(/chem_surfvals_get('N2OVMR')/))
            ierr=pio_put_var (tape(t)%File, tape(t)%f11vmrid,(/start/), (/count1/),(/chem_surfvals_get('F11VMR')/))
            ierr=pio_put_var (tape(t)%File, tape(t)%f12vmrid,(/start/), (/count1/),(/chem_surfvals_get('F12VMR')/))
            ierr=pio_put_var (tape(t)%File, tape(t)%sol_tsiid,(/start/), (/count1/),(/sol_tsi/))

            if (solar_parms_on) then
              ierr=pio_put_var (tape(t)%File, tape(t)%f107id, (/start/), (/count1/),(/ f107 /) )
              ierr=pio_put_var (tape(t)%File, tape(t)%f107aid,(/start/), (/count1/),(/ f107a /) )
              ierr=pio_put_var (tape(t)%File, tape(t)%f107pid,(/start/), (/count1/),(/ f107p /) )
              ierr=pio_put_var (tape(t)%File, tape(t)%kpid,   (/start/), (/count1/),(/ kp /) )
              ierr=pio_put_var (tape(t)%File, tape(t)%apid,   (/start/), (/count1/),(/ ap /) )
            endif
            if (solar_wind_on) then
              ierr=pio_put_var (tape(t)%File, tape(t)%byimfid, (/start/), (/count1/),(/ byimf /) )
              ierr=pio_put_var (tape(t)%File, tape(t)%bzimfid, (/start/), (/count1/),(/ bzimf /) )
              ierr=pio_put_var (tape(t)%File, tape(t)%swvelid, (/start/), (/count1/),(/ swvel /) )
              ierr=pio_put_var (tape(t)%File, tape(t)%swdenid, (/start/), (/count1/),(/ swden /) )
            endif
          end if

          ierr = pio_put_var (tape(t)%File, tape(t)%datesecid,(/start/),(/count1/),(/ncsec/))
#if ( defined BFB_CAM_SCAM_IOP )
          dtime = get_step_size()
          tsec=dtime*nstep
          ierr = pio_put_var (tape(t)%File, tape(t)%tsecid,(/start/),(/count1/),(/tsec/))
#endif
          ierr = pio_put_var (tape(t)%File, tape(t)%nstephid,(/start/),(/count1/),(/nstep/))
          time = ndcur + nscur/86400._r8
          ierr=pio_put_var (tape(t)%File, tape(t)%timeid, (/start/),(/count1/),(/time/))

          startc(1) = 1
          startc(2) = start
          countc(1) = 2
          countc(2) = 1
          if (is_initfile(file_index=t)) then
            tdata = time   ! Inithist file is always instantanious data
          else
            tdata(1) = beg_time(t)
            tdata(2) = time
          end if
          ierr=pio_put_var (tape(t)%File, tape(t)%tbndid, startc, countc, tdata)
          if(.not.restart) beg_time(t) = time  ! update beginning time of next interval
          startc(1) = 1
          startc(2) = start
          countc(1) = 8
          countc(2) = 1
          call datetime (cdate, ctime)
          ierr = pio_put_var (tape(t)%File, tape(t)%date_writtenid, startc, countc, (/cdate/))
          ierr = pio_put_var (tape(t)%File, tape(t)%time_writtenid, startc, countc, (/ctime/))

          if(.not. restart) then
            !$OMP PARALLEL DO PRIVATE (F)
            do f=1,nflds(t)
              ! Normalized averaged fields
              if (tape(t)%hlist(f)%avgflag /= 'I') then
                call h_normalize (f, t)
              end if
            end do
          end if
          !
          ! Write field to history tape.  Note that this is NOT threaded due to netcdf limitations
          !
          call t_startf ('dump_field')
          do f=1,nflds(t)
            call dump_field(f, t, restart)
          end do
          call t_stopf ('dump_field')
          !
          ! Zero history buffers and accumulators now that the fields have been written.
          !



          if(restart) then
            do f=1,nflds(t)
              if(associated(tape(t)%hlist(f)%varid)) then
                deallocate(tape(t)%hlist(f)%varid)
                nullify(tape(t)%hlist(f)%varid)
              end if
            end do
            call cam_pio_closefile(tape(t)%File)
          else
            !$OMP PARALLEL DO PRIVATE (F)
            do f=1,nflds(t)
              call h_zero (f, t)
            end do
          end if
        end if
      end if

    end do

    return
  end subroutine wshist

  !#######################################################################

  subroutine addfld_1d(fname, vdim_name, avgflag, units, long_name,           &
       gridname, flag_xyfill, sampling_seq, standard_name, fill_value)

    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Add a field to the master field list
    !
    ! Method: Put input arguments of field name, units, number of levels,
    !         averaging flag, and long name into a type entry in the global
    !         master field list (masterlist).
    !
    !-----------------------------------------------------------------------

    !
    ! Arguments
    !
    character(len=*), intent(in)  :: fname      ! field name (max_fieldname_len)
    character(len=*), intent(in)  :: vdim_name  ! NetCDF dimension name (or scalar coordinate)
    character(len=1), intent(in)  :: avgflag    ! averaging flag
    character(len=*), intent(in)  :: units      ! units of fname (max_chars)
    character(len=*), intent(in)  :: long_name  ! long name of field (max_chars)

    character(len=*), intent(in), optional :: gridname    ! decomposition type
    logical, intent(in), optional :: flag_xyfill ! non-applicable xy points flagged with fillvalue
    character(len=*), intent(in), optional :: sampling_seq ! sampling sequence - if not every timestep,
    ! how often field is sampled:
    ! every other; only during LW/SW radiation calcs, etc.
    character(len=*), intent(in), optional :: standard_name  ! CF standard name (max_chars)
    real(r8),         intent(in), optional :: fill_value

    !
    ! Local workspace
    !
    character(len=max_chars), allocatable :: dimnames(:)
    integer                               :: index

    if (trim(vdim_name) == trim(horiz_only)) then
      allocate(dimnames(0))
    else
      index = get_hist_coord_index(trim(vdim_name))
      if (index < 1) then
        call endrun('ADDFLD: Invalid coordinate, '//trim(vdim_name))
      end if
      allocate(dimnames(1))
      dimnames(1) = trim(vdim_name)
    end if
    call addfld(fname, dimnames, avgflag, units, long_name, gridname,         &
         flag_xyfill, sampling_seq, standard_name, fill_value)

  end subroutine addfld_1d

  subroutine addfld_nd(fname, dimnames, avgflag, units, long_name,            &
       gridname, flag_xyfill, sampling_seq, standard_name, fill_value)

    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Add a field to the master field list
    !
    ! Method: Put input arguments of field name, units, number of levels,
    !         averaging flag, and long name into a type entry in the global
    !         master field list (masterlist).
    !
    !-----------------------------------------------------------------------
    use cam_history_support, only: fillvalue, hist_coord_find_levels
    use cam_grid_support,    only: cam_grid_id, cam_grid_is_zonal
    use cam_grid_support,    only: cam_grid_get_coord_names
    use constituents,        only: pcnst, cnst_get_ind, cnst_get_type_byind

    !
    ! Arguments
    !
    character(len=*), intent(in)  :: fname      ! field name (max_fieldname_len)
    character(len=*), intent(in)  :: dimnames(:) ! NetCDF dimension names (except grid dims)
    character(len=1), intent(in)  :: avgflag    ! averaging flag
    character(len=*), intent(in)  :: units      ! units of fname (max_chars)
    character(len=*), intent(in)  :: long_name  ! long name of field (max_chars)

    character(len=*), intent(in), optional :: gridname    ! decomposition type
    logical, intent(in), optional :: flag_xyfill ! non-applicable xy points flagged with fillvalue
    character(len=*), intent(in), optional :: sampling_seq ! sampling sequence - if not every timestep,
    ! how often field is sampled:
    ! every other; only during LW/SW radiation calcs, etc.
    character(len=*), intent(in), optional :: standard_name  ! CF standard name (max_chars)
    real(r8),         intent(in), optional :: fill_value

    !
    ! Local workspace
    !
    character(len=max_fieldname_len) :: fname_tmp ! local copy of fname
    character(len=max_fieldname_len) :: coord_name ! for cell_methods
    character(len=128)               :: errormsg
    character(len=3)                 :: mixing_ratio
    type(master_entry), pointer      :: listentry

    integer :: dimcnt
    integer :: idx

    if (htapes_defined) then
      call endrun ('ADDFLD: Attempt to add field '//trim(fname)//' after history files set')
    end if

    !
    ! Ensure that new field name is not all blanks
    !
    if (len_trim(fname)==0) then
      call endrun('ADDFLD: blank field name not allowed')
    end if
    !
    ! Ensure that new field name is not longer than allowed
    ! (strip "&IC" suffix if it exists)
    !
    fname_tmp  = fname
    fname_tmp  = strip_suffix(fname_tmp)

    if (len_trim(fname_tmp) > fieldname_len) then
      write(iulog,*)'ADDFLD: field name cannot be longer than ',fieldname_len,' characters long'
      write(iulog,*)'Field name:  ',fname
      write(errormsg, *) 'Field name, "', trim(fname), '" is too long'
      call endrun('ADDFLD: '//trim(errormsg))
    end if
    !
    ! Ensure that new field doesn't already exist
    !
    listentry => get_entry_by_name(masterlinkedlist, fname)
    if(associated(listentry)) then
      call endrun ('ADDFLD:  '//fname//' already on list')
    end if

    ! If the field is an advected constituent determine whether its concentration
    ! is based on dry or wet air.
    call cnst_get_ind(fname_tmp, idx, abort=.false.)
    mixing_ratio = ''
    if (idx > 0) then
       mixing_ratio = cnst_get_type_byind(idx)
    end if

    ! Add field to Master Field List arrays fieldn and iflds
    !
    allocate(listentry)
    listentry%field%name         = fname
    listentry%field%long_name    = long_name
    listentry%field%numlev       = 1        ! Will change if lev or ilev in shape
    listentry%field%units        = units
    listentry%field%mixing_ratio = mixing_ratio
    listentry%field%meridional_complement = -1
    listentry%field%zonal_complement      = -1
    listentry%htapeindx(:) = -1
    listentry%act_sometape = .false.
    listentry%actflag(:) = .false.

    ! Make sure we have a valid gridname
    if (present(gridname)) then
      listentry%field%decomp_type = cam_grid_id(trim(gridname))
    else
      listentry%field%decomp_type = cam_grid_id('physgrid')
    end if
    if (listentry%field%decomp_type < 0) then
      write(errormsg, *) 'Invalid grid name, "', trim(gridname), '" for ',    &
           trim(fname)
      call endrun('ADDFLD: '//trim(errormsg))
    end if

    !
    ! Indicate sampling sequence of field (i.e., how often "outfld" is called)
    ! If not every timestep (default), then give a descriptor indicating the
    ! sampling pattern.  Currently, the only valid value is "rad_lwsw" for sampling
    ! during LW/SW radiation timesteps only
    !
    if (present(sampling_seq)) then
      listentry%field%sampling_seq = sampling_seq
    else
      listentry%field%sampling_seq = ' '
    end if
    ! Indicate if some field pre-processing occurred (e.g., zonal mean)
    if (cam_grid_is_zonal(listentry%field%decomp_type)) then
      call cam_grid_get_coord_names(listentry%field%decomp_type, coord_name, errormsg)
      ! Zonal method currently hardcoded to 'mean'.
      listentry%field%cell_methods = trim(coord_name)//': mean'
    else
      listentry%field%cell_methods = ''
    end if
    !
    ! Whether to apply xy fillvalue: default is false
    !
    if (present(flag_xyfill)) then
      listentry%field%flag_xyfill = flag_xyfill
    else
      listentry%field%flag_xyfill = .false.
    end if

    !
    !    Allow external packages to have fillvalues different than default
    !

    if(present(fill_value)) then
      listentry%field%fillvalue = fill_value
    else
      listentry%field%fillvalue = fillvalue
    endif

    !
    ! Process shape
    !

    if (associated(listentry%field%mdims)) then
      deallocate(listentry%field%mdims)
    end if
    nullify(listentry%field%mdims)
    dimcnt = size(dimnames)
    allocate(listentry%field%mdims(dimcnt))
    call lookup_hist_coord_indices(dimnames, listentry%field%mdims)
    if(dimcnt > maxvarmdims) then
      maxvarmdims = dimcnt
    end if
    ! Check for subcols (currently limited to first dimension)
    listentry%field%is_subcol = .false.
    if (size(dimnames) > 0) then
      if (trim(dimnames(1)) == 'psubcols') then
        if (listentry%field%decomp_type /= cam_grid_id('physgrid')) then
          write(errormsg, *) "Cannot add ", trim(fname),                      &
               "Subcolumn history output only allowed on physgrid"
          call endrun("ADDFLD: "//errormsg)
          listentry%field%is_subcol = .true.
        end if
      end if
    end if
    ! Levels
    listentry%field%numlev = hist_coord_find_levels(dimnames)
    if (listentry%field%numlev <= 0) then
      listentry%field%numlev = 1
    end if

    !
    ! Dimension history info based on decomposition type (grid)
    !
    call set_field_dimensions(listentry%field)

    !
    ! These 2 fields are used only in master field list, not runtime field list
    !
    listentry%avgflag(:) = avgflag
    listentry%actflag(:) = .false.

    do dimcnt = 1, ptapes
      call AvgflagToString(avgflag, listentry%time_op(dimcnt))
    end do

    nullify(listentry%next_entry)

    call add_entry_to_master(listentry)
    return
  end subroutine addfld_nd

  !#######################################################################

  ! field_part_of_vector: Determinie if fname is part of a vector set
  !       Optionally fill in the names of the vector set fields
  logical function field_part_of_vector(fname, meridional_name, zonal_name)

    ! Dummy arguments
    character(len=*),           intent(in)  :: fname
    character(len=*), optional, intent(out) :: meridional_name
    character(len=*), optional, intent(out) :: zonal_name

    ! Local variables
    type(master_entry), pointer             :: listentry

    listentry => get_entry_by_name(masterlinkedlist, fname)
    if (associated(listentry)) then
      if ( (len_trim(listentry%meridional_field) > 0) .or.                     &
           (len_trim(listentry%zonal_field) > 0)) then
        field_part_of_vector = .true.
        if (present(meridional_name)) then
          meridional_name = listentry%meridional_field
        end if
        if (present(zonal_name)) then
          zonal_name = listentry%zonal_field
        end if
      else
        field_part_of_vector = .false.
      end if
    else
      field_part_of_vector = .false.
    end if
    if (.not. field_part_of_vector) then
      if (present(meridional_name)) then
        meridional_name = ''
      end if
      if (present(zonal_name)) then
        zonal_name = ''
      end if
    end if

  end function field_part_of_vector


  ! register_vector_field: Register a pair of history field names as
  !           being a vector complement set.
  !           This information is used to set up interpolated history output.
  ! NB: register_vector_field must be called after both fields are defined
  !     with addfld
  subroutine register_vector_field(zonal_field_name, meridional_field_name)

    ! Dummy arguments
    character(len=*),             intent(in) :: zonal_field_name
    character(len=*),             intent(in) :: meridional_field_name

    ! Local variables
    type(master_entry), pointer      :: mlistentry
    type(master_entry), pointer      :: zlistentry
    character(len=*),   parameter    :: subname = 'REGISTER_VECTOR_FIELD'
    character(len=max_chars)         :: errormsg

    if (htapes_defined) then
      write(errormsg, '(5a)') ': Attempt to register vector field (',         &
           trim(zonal_field_name), ', ', trim(meridional_field_name),         &
           ') after history files set'
      call endrun (trim(subname)//errormsg)
    end if

    ! Look for the field IDs
    zlistentry => get_entry_by_name(masterlinkedlist, zonal_field_name)
    mlistentry => get_entry_by_name(masterlinkedlist, meridional_field_name)
    ! Has either of these fields been previously registered?
    if (associated(mlistentry)) then
      if (len_trim(mlistentry%meridional_field) > 0) then
        write(errormsg, '(9a)') ': ERROR attempting to register vector ',     &
             'field (', trim(zonal_field_name), ', ',                         &
             trim(meridional_field_name), '), ', trim(meridional_field_name), &
             ' has been registered as part of a vector field with ',          &
             trim(mlistentry%meridional_field)
        call endrun (trim(subname)//errormsg)
      else if (len_trim(mlistentry%zonal_field) > 0) then
        write(errormsg, '(9a)') ': ERROR attempting to register vector ',     &
             'field (', trim(zonal_field_name), ', ',                         &
             trim(meridional_field_name), '), ', trim(meridional_field_name), &
             ' has been registered as part of a vector field with ',          &
             trim(mlistentry%zonal_field)
        call endrun (trim(subname)//errormsg)
      end if
    end if
    if (associated(zlistentry)) then
      if (len_trim(zlistentry%meridional_field) > 0) then
        write(errormsg, '(9a)') ': ERROR attempting to register vector ',     &
             'field (', trim(zonal_field_name), ', ',                         &
             trim(meridional_field_name), '), ', trim(zonal_field_name),      &
             ' has been registered as part of a vector field with ',          &
             trim(zlistentry%meridional_field)
        call endrun (trim(subname)//errormsg)
      else if (len_trim(zlistentry%zonal_field) > 0) then
        write(errormsg, '(9a)') ': ERROR attempting to register vector ',     &
             'field (', trim(zonal_field_name), ', ',                         &
             trim(meridional_field_name), '), ', trim(zonal_field_name),      &
             ' has been registered as part of a vector field with ',          &
             trim(zlistentry%meridional_field)
        call endrun (trim(subname)//errormsg)
      end if
    end if
    if(associated(mlistentry) .and. associated(zlistentry)) then
      zlistentry%meridional_field = mlistentry%field%name
      zlistentry%zonal_field      = ''
      mlistentry%meridional_field = ''
      mlistentry%zonal_field      = zlistentry%field%name
    else if (associated(mlistentry)) then
      write(errormsg, '(7a)') ': ERROR attempting to register vector field (',&
           trim(zonal_field_name), ', ', trim(meridional_field_name),         &
           '), ', trim(zonal_field_name), ' is not defined'
      call endrun (trim(subname)//errormsg)
    else if (associated(zlistentry)) then
      write(errormsg, '(7a)') ': ERROR attempting to register vector field (',&
           trim(zonal_field_name), ', ', trim(meridional_field_name),         &
           '), ', trim(meridional_field_name), ' is not defined'
      call endrun (trim(subname)//errormsg)
    else
      write(errormsg, '(5a)') ': ERROR attempting to register vector field (',&
           trim(zonal_field_name), ', ', trim(meridional_field_name),         &
           '), neither field is defined'
      call endrun (trim(subname)//errormsg)
    end if
  end subroutine register_vector_field

  subroutine add_entry_to_master( newentry)
    type(master_entry), target, intent(in) :: newentry
    type(master_entry), pointer :: listentry

    if(associated(masterlinkedlist)) then
      listentry => masterlinkedlist
      do while(associated(listentry%next_entry))
        listentry=>listentry%next_entry
      end do
      listentry%next_entry=>newentry
    else
      masterlinkedlist=>newentry
    end if

  end subroutine add_entry_to_master

  !#######################################################################

  subroutine wrapup (rstwr, nlend)
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose:
    ! Close history files.
    !
    ! Method:
    ! This routine will close any full hist. files
    ! or any hist. file that has data on it when restart files are being
    ! written.
    ! If a partially full history file was disposed (for restart
    ! purposes), then wrapup will open that unit back up and position
    ! it for appending new data.
    !
    ! Original version: CCM2
    !
    !-----------------------------------------------------------------------
    !
    use pio, only : pio_file_is_open
    use shr_kind_mod,  only: r8 => shr_kind_r8
    use ioFileMod
    use time_manager,  only: get_nstep, get_curr_date, get_curr_time
    use cam_pio_utils, only: cam_pio_openfile, cam_pio_closefile

    !
    ! Input arguments
    !
    logical, intent(in) :: rstwr   ! true => restart files are written this timestep
    logical, intent(in) :: nlend   ! Flag if time to end

    !
    ! Local workspace
    !
    integer  :: nstep            ! current timestep number
    integer  :: ncsec            ! time of day relative to current date [secs]
    integer  :: ndcur            ! days component of current time
    integer  :: nscur            ! seconds component of current time
    integer  :: yr, mon, day     ! year, month, day components of a date

    logical  :: lfill   (ptapes) ! Is history file ready to dispose?
    logical  :: lhdisp           ! true => history file is disposed
    logical  :: lhfill           ! true => history file is full

    integer  :: t                ! History file number
    integer  :: f
    real(r8) :: tday             ! Model day number for printout
    !-----------------------------------------------------------------------

    tape => history_tape

    nstep = get_nstep()
    call get_curr_date(yr, mon, day, ncsec)
    call get_curr_time(ndcur, nscur)
    !
    !-----------------------------------------------------------------------
    ! Dispose history files.
    !-----------------------------------------------------------------------
    !
    ! Begin loop over ptapes (the no. of declared history files - primary
    ! and auxiliary).  This loop disposes a history file to Mass Store
    ! when appropriate.
    !
    do t=1,ptapes
      if (nflds(t) == 0) cycle

      lfill(t) = .false.
      !
      ! Find out if file is full
      !
      if (hstwr(t) .and. nfils(t) >= mfilt(t)) then
        lfill(t) = .true.
      endif
      !
      ! Dispose history file if
      !    1) file is filled or
      !    2) this is the end of run and file has data on it or
      !    3) restarts are being put out and history file has data on it
      !
      if (lfill(t) .or. (nlend .and. nfils(t) >= 1) .or. (rstwr .and. nfils(t) >= 1)) then
        !
        ! Dispose history file
        !
        !
        ! Is this the 0 timestep data of a monthly run?
        ! If so, just close primary unit do not dispose.
        !
        if (masterproc) write(iulog,*)'WRAPUP: nf_close(',t,')=',trim(nhfil(t))
        if(pio_file_is_open(tape(t)%File)) then
          if (nlend .or. lfill(t)) then
            do f=1,nflds(t)
              if (associated(tape(t)%hlist(f)%varid)) then
                deallocate(tape(t)%hlist(f)%varid)
                nullify(tape(t)%hlist(f)%varid)
              end if
            end do
          end if
          call cam_pio_closefile(tape(t)%File)
        end if
        if (nhtfrq(t) /= 0 .or. nstep > 0) then

          !
          ! Print information concerning model output.
          ! Model day number = iteration number of history file data * delta-t / (seconds per day)
          !
          tday = ndcur + nscur/86400._r8
          if(masterproc) then
            if (t==1) then
              write(iulog,*)'   Primary history file'
            else
              write(iulog,*)'   Auxiliary history file number ', t-1
            end if
            write(iulog,9003)nstep,nfils(t),tday
            write(iulog,9004)
          end if
          !
          ! Auxilary files may have been closed and saved off without being full.
          ! We must reopen the files and position them for more data.
          ! Must position auxiliary files if not full
          !
          if (.not.nlend .and. .not.lfill(t)) then
            call cam_PIO_openfile (tape(t)%File, nhfil(t), PIO_WRITE)
            call h_inquire(t)
          end if
        endif                 ! if 0 timestep of montly run****
      end if                      ! if time dispose history fiels***
    end do                         ! do ptapes
    !
    ! Reset number of files on each history tape
    !
    do t=1,ptapes
      if (nflds(t) == 0) cycle
      lhfill = hstwr(t) .and. nfils(t) >= mfilt(t)
      lhdisp = lhfill .or. (nlend .and. nfils(t) >= 1) .or. &
           (rstwr .and. nfils(t) >= 1)
      if (lhfill.and.lhdisp) then
        nfils(t) = 0
      endif
    end do
    return
9003 format('    Output at NSTEP     = ',i10,/, &
         '    Number of time samples on this file = ',i10,/, &
         '    Model Day           = ',f10.2)
9004 format('---------------------------------------')
  end subroutine wrapup


  integer function gen_hash_key(string)
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Generate a hash key on the interval [0 .. tbl_hash_pri_sz-1]
    !          given a character string.
    !
    ! Algorithm is a variant of perl's internal hashing function.
    !
    !-----------------------------------------------------------------------
    !
    implicit none
    !
    !  Arguments:
    !
    character(len=*), intent(in) :: string
    !
    !  Local.
    !
    integer :: hash
    integer :: i

    hash = gen_hash_key_offset

    if ( len(string) /= 19 ) then
      !
      !     Process arbitrary string length.
      !
      do i = 1, len(string)
        hash = ieor(hash, (ichar(string(i:i)) * tbl_gen_hash_key(iand(i-1,tbl_max_idx))))
      end do
    else
      !
      !     Special case string length = 19
      !
      hash = ieor(hash , ichar(string(1:1))   * 61)
      hash = ieor(hash , ichar(string(2:2))   * 59)
      hash = ieor(hash , ichar(string(3:3))   * 53)
      hash = ieor(hash , ichar(string(4:4))   * 47)
      hash = ieor(hash , ichar(string(5:5))   * 43)
      hash = ieor(hash , ichar(string(6:6))   * 41)
      hash = ieor(hash , ichar(string(7:7))   * 37)
      hash = ieor(hash , ichar(string(8:8))   * 31)
      hash = ieor(hash , ichar(string(9:9))   * 29)
      hash = ieor(hash , ichar(string(10:10)) * 23)
      hash = ieor(hash , ichar(string(11:11)) * 17)
      hash = ieor(hash , ichar(string(12:12)) * 13)
      hash = ieor(hash , ichar(string(13:13)) * 11)
      hash = ieor(hash , ichar(string(14:14)) * 7)
      hash = ieor(hash , ichar(string(15:15)) * 3)
      hash = ieor(hash , ichar(string(16:16)) * 1)
      hash = ieor(hash , ichar(string(17:17)) * 61)
      hash = ieor(hash , ichar(string(18:18)) * 59)
      hash = ieor(hash , ichar(string(19:19)) * 53)
    end if

    gen_hash_key = iand(hash, tbl_hash_pri_sz-1)

    return

  end function gen_hash_key

  !#######################################################################

  integer function get_masterlist_indx(fldname)
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Return the the index of the field's name on the master file list.
    !
    !          If the field is not found on the masterlist, return -1.
    !
    !-----------------------------------------------------------------------
    !
    !  Arguments:
    !
    character(len=*), intent(in) :: fldname
    !
    !  Local.
    !
    integer :: hash_key
    integer :: ff
    integer :: ii
    integer :: io   ! Index of overflow chain in overflow table
    integer :: in   ! Number of entries on overflow chain

    hash_key = gen_hash_key(fldname)
    ff = tbl_hash_pri(hash_key)
    if ( ff < 0 ) then
      io = abs(ff)
      in = tbl_hash_oflow(io)
      do ii = 1, in
        ff = tbl_hash_oflow(io+ii)
        if ( masterlist(ff)%thisentry%field%name == fldname ) exit
      end do
    end if

    if (ff == 0) then
      ! fldname generated a hash key that doesn't have an entry in tbl_hash_pri.
      ! This means that fldname isn't in the masterlist
      get_masterlist_indx = -1
      return
      call endrun ('GET_MASTERLIST_INDX: attemping to output field '//fldname//' not on master list')
    end if

    if (associated(masterlist(ff)%thisentry) .and. masterlist(ff)%thisentry%field%name /= fldname ) then
      get_masterlist_indx = -1
      return
      call endrun ('GET_MASTERLIST_INDX: error finding field '//fldname//' on master list')
    end if

    get_masterlist_indx = ff
    return
  end function get_masterlist_indx
  !#######################################################################

  subroutine bld_outfld_hash_tbls()
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Build primary and overflow hash tables for outfld processing.
    !
    ! Steps:
    !  1) Foreach field on masterlist, find all collisions.
    !  2) Given the number of collisions, verify overflow table has sufficient
    !     space.
    !  3) Build primary and overflow indices.
    !
    !-----------------------------------------------------------------------
    !
    !  Local.
    !
    integer :: ff
    integer :: ii
    integer :: itemp
    integer :: ncollisions
    integer :: hash_key
    type(master_entry), pointer :: listentry
    !
    !  1) Find all collisions.
    !
    tbl_hash_pri = 0

    ff=0
    allocate(masterlist(nfmaster))
    listentry=>masterlinkedlist
    do while(associated(listentry))
      ff=ff+1
      masterlist(ff)%thisentry=>listentry
      listentry=>listentry%next_entry
    end do
    if(ff /= nfmaster) then
      write(iulog,*) 'nfmaster = ',nfmaster, ' ff=',ff
      call endrun('mismatch in expected size of nfmaster')
    end if


    do ff = 1, nfmaster
      hash_key = gen_hash_key(masterlist(ff)%thisentry%field%name)
      tbl_hash_pri(hash_key) = tbl_hash_pri(hash_key) + 1
    end do

    !
    !  2) Count number of collisions and define start of a individual
    !     collision's chain in overflow table. A collision is defined to be any
    !     location in tbl_hash_pri that has a value > 1.
    !
    ncollisions = 0
    do ii = 0, tbl_hash_pri_sz-1
      if ( tbl_hash_pri(ii) > 1 ) then  ! Define start of chain in O.F. table
        itemp = tbl_hash_pri(ii)
        tbl_hash_pri(ii) = -(ncollisions + 1)
        ncollisions = ncollisions + itemp + 1
      end if
    end do

    if ( ncollisions > tbl_hash_oflow_sz ) then
      write(iulog,*) 'BLD_OUTFLD_HASH_TBLS: ncollisions > tbl_hash_oflow_sz', &
           ncollisions, tbl_hash_oflow_sz
      call endrun()
    end if

    !
    !  3) Build primary and overflow tables.
    !     i - set collisions in tbl_hash_pri to point to their respective
    !         chain in the overflow table.
    !
    tbl_hash_oflow = 0

    do ff = 1, nfmaster
      hash_key = gen_hash_key(masterlist(ff)%thisentry%field%name)
      if ( tbl_hash_pri(hash_key) < 0 ) then
        ii = abs(tbl_hash_pri(hash_key))
        tbl_hash_oflow(ii) = tbl_hash_oflow(ii) + 1
        tbl_hash_oflow(ii+tbl_hash_oflow(ii)) = ff
      else
        tbl_hash_pri(hash_key) = ff
      end if
    end do

    !
    !  Dump out primary and overflow hashing tables.
    !
    !   if ( masterproc ) then
    !      do ii = 0, tbl_hash_pri_sz-1
    !         if ( tbl_hash_pri(ii) /= 0 ) write(iulog,666) 'tbl_hash_pri', ii, tbl_hash_pri(ii)
    !      end do
    !
    !      do ii = 1, tbl_hash_oflow_sz
    !         if ( tbl_hash_oflow(ii) /= 0 ) write(iulog,666) 'tbl_hash_oflow', ii, tbl_hash_oflow(ii)
    !      end do
    !
    !      itemp = 0
    !      ii = 1
    !      do
    !         if ( tbl_hash_oflow(ii) == 0 ) exit
    !         itemp = itemp + 1
    !         write(iulog,*) 'Overflow chain ', itemp, ' has ', tbl_hash_oflow(ii), ' entries:'
    !         do ff = 1, tbl_hash_oflow(ii)  ! dump out colliding names on this chain
    !            write(iulog,*) '     ', ff, ' = ', tbl_hash_oflow(ii+ff), &
    !                       ' ', masterlist(tbl_hash_oflow(ii+ff))%thisentry%field%name
    !         end do
    !         ii = ii + tbl_hash_oflow(ii) +1 !advance pointer to start of next chain
    !      end do
    !   end if

    return
666 format(1x, a, '(', i4, ')', 1x, i6)

  end subroutine bld_outfld_hash_tbls

  !#######################################################################

  subroutine bld_htapefld_indices
    !
    !-----------------------------------------------------------------------
    !
    ! Purpose: Set history tape field indicies in masterlist for each
    !          field defined on every tape.
    !
    ! Note: because of restart processing, the actflag field is cleared and
    !       then set only for active output fields on the different history
    !       tapes.
    !
    !-----------------------------------------------------------------------
    !
    !  Arguments:
    !

    !
    !  Local.
    !
    integer :: f
    integer :: t

    !
    !  Initialize htapeindx to an invalid value.
    !
    type(master_entry), pointer :: listentry

    ! reset all the active flags to false
    ! this is needed so that restarts work properly -- fvitt
    listentry=>masterlinkedlist
    do while(associated(listentry))
      listentry%actflag(:) = .false.
      listentry%act_sometape = .false.
      listentry=>listentry%next_entry
    end do

    do t = 1, ptapes
      do f = 1, nflds(t)
        listentry => get_entry_by_name(masterlinkedlist, tape(t)%hlist(f)%field%name)
        if(.not.associated(listentry)) then
          write(iulog,*) 'BLD_HTAPEFLD_INDICES: something wrong, field not found on masterlist'
          write(iulog,*) 'BLD_HTAPEFLD_INDICES: t, f, ff = ', t, f
          write(iulog,*) 'BLD_HTAPEFLD_INDICES: tape%name = ', tape(t)%hlist(f)%field%name
          call endrun
        end if
        listentry%act_sometape = .true.
        listentry%actflag(t) = .true.
        listentry%htapeindx(t) = f
      end do
    end do

    !
    ! set flag indicating h-tape contents are now defined (needed by addfld)
    !
    htapes_defined = .true.

    return
  end subroutine bld_htapefld_indices

  !#######################################################################

  logical function hist_fld_active(fname)
    !
    !------------------------------------------------------------------------
    !
    ! Purpose: determine if a field is active on any history file
    !
    !------------------------------------------------------------------------
    !
    ! Arguments
    !
    character(len=*), intent(in) :: fname ! Field name
    !
    ! Local variables
    !
    character*(max_fieldname_len) :: fname_loc  ! max-char equivalent of fname
    integer :: ff                  ! masterlist index pointer
    !-----------------------------------------------------------------------

    fname_loc = fname
    ff = get_masterlist_indx(fname_loc)
    if ( ff < 0 ) then
      hist_fld_active = .false.
    else
      hist_fld_active = masterlist(ff)%thisentry%act_sometape
    end if

  end function hist_fld_active

  !#######################################################################

  function hist_fld_col_active(fname, lchnk, numcols)
    use cam_history_support, only: history_patch_t

    ! Determine whether each column in a field is active on any history file.
    ! The purpose of this routine is to provide information which would allow
    ! a diagnostic physics parameterization to only be run on a subset of
    ! columns in the case when only column or regional output is requested.
    !
    ! **N.B.** The field is assumed to be using the physics decomposition.

    ! Arguments
    character(len=*), intent(in) :: fname   ! Field name
    integer,          intent(in) :: lchnk   ! chunk ID
    integer,          intent(in) :: numcols ! Size of return array

    ! Return value
    logical :: hist_fld_col_active(numcols)

    ! Local variables
    integer                         :: ff          ! masterlist index pointer
    integer                         :: i
    integer                         :: t           ! history file (tape) index
    integer                         :: f           ! field index
    integer                         :: decomp
    logical                         :: activeloc(numcols)
    integer                         :: num_patches
    logical                         :: patch_output
    logical                         :: found
    type(history_patch_t), pointer  :: patchptr

    type (active_entry),   pointer  :: tape(:)

    !-----------------------------------------------------------------------

    ! Initialize to false.  Then look to see if and where active.
    hist_fld_col_active = .false.

    ! Check for name in the master list.
    call get_field_properties(fname, found, tape_out=tape, ff_out=ff)

    ! If not in master list then return.
    if (.not. found) return

    ! If in master list, but not active on any file then return
    if (.not. masterlist(ff)%thisentry%act_sometape) return

    ! Loop over history files and check for the field/column in each one
    do t = 1, ptapes

      ! Is the field active in this file?  If not the cycle to next file.
      if (.not. masterlist(ff)%thisentry%actflag(t)) cycle

      f = masterlist(ff)%thisentry%htapeindx(t)
      decomp = tape(t)%hlist(f)%field%decomp_type
      patch_output = associated(tape(t)%patches)

      ! Check whether this file has patch (column) output.
      if (patch_output) then
        num_patches = size(tape(t)%patches)

        do i = 1, num_patches
          patchptr => tape(t)%patches(i)
          activeloc = .false.
          call patchptr%active_cols(decomp, lchnk, activeloc)
          hist_fld_col_active = hist_fld_col_active .or. activeloc
        end do
      else

        ! No column output has been requested.  In that case the field has
        ! global output which implies all columns are active.  No need to
        ! check any other history files.
        hist_fld_col_active = .true.
        exit

      end if

    end do ! history files

  end function hist_fld_col_active

  subroutine cam_history_snapshot_deactivate(name)

  ! This subroutine deactivates (sets actflag to false) for all tapes

  character(len=*), intent(in) :: name

  logical :: found
  integer :: ff

  call get_field_properties(trim(name), found, ff_out=ff, no_tape_check_in=.true.)
  masterlist(ff)%thisentry%actflag(:) = .false.

  end subroutine cam_history_snapshot_deactivate

  subroutine cam_history_snapshot_activate(name, tape)

  ! This subroutine activates (set aftflag to true) for the requested tape number

  character(len=*), intent(in) :: name
  integer,          intent(in) :: tape

  logical :: found
  integer :: ff

  call get_field_properties(trim(name), found, ff_out=ff, no_tape_check_in=.true.)
  masterlist(ff)%thisentry%actflag(tape) = .true.

  end subroutine cam_history_snapshot_activate

end module cam_history
